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"  A 

This  dissertation  addresses  the  problem  of  determining 
the  correct  relationship  between  the  statistics  of  a  con¬ 
tinuous  random  process  and  the  statistics  of  a  discrete 
random  process  used  to  simulate  the  continuous  random  pro¬ 
cess.  The  findings  of  this  research  are  directly  applica¬ 
ble  to  the  general  field  of  digital  simulation  of  physical 
systems  described  by  ordinary  differential  equations. 

It  is  shown  that  to  ensure  a  faithful  digital  simula¬ 
tion  of  a  continuous  random  process,  the  noise  statistics 

V» 

of  the  random  number  generator  must  be  set  to  values  dras¬ 
tically  different  from  the  noise  statistics  of  the  contin¬ 
uous  random  process.  Further,  it  is  established  that  the 


relationship  between  the  continuous  and  discrete  statistics 
is  a  function  of  the  integration  method  used  in  the  digital 
simulation. 

The  proper  functional  relationship  between  the  dis¬ 
crete  and  continuous  noise  statistics  is  derived  for 

1 —  the  class  of  Runge-Kutta  integrators, 

2.  the  4th  order  Adams-Bashforth  integrator,  and 

3.  the  Adams-Moul ton  corrector  formula. 

The  derived  relationships  are  applied  to  a  specific  problem 
and  are  demonstrated  by  simulation.  The  simulation  results 
are  compared  to  exact  solutions.  Additionally,  the  require¬ 
ment  for  proper  operation  of  a  variable-step-size  algorithm 
is  developed. 
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I.  INTRODUCTION 


The  design  of  modern  technological  systems  is  a  well 
defined  logical  process  involving  many  phases  beginning 
with  a  statement  of  need  and  culminating  in  hardware.  The 
success  of  the  final  product  is  critically  dependent  on  the 
successful  completion  of  each  design  phase.  In  each  of  the 
phases,  analysis  plays  a  fundamental  role,  whether  the  goal 
of  the  analysis  is  to  test  the  potential  payoff  of  a  feasi¬ 
ble  solution,  or  simply  to  further  the  engineer's  under¬ 
standing  of  the  problem.  The  importance  of  an  accurate 
analysis,  regardless  of  the  design  phase  for  which  it  is 
performed,  can  not  be  overemphasized. 

To  perform  an  accurate  analysis,  the  analyst  has  a 
number  of  useful  and  time-proven  tools.  One  of  the  most 
useful  tools  for  analysis  is  based  on  the  modern  engineer's 
ability  to  describe  a  physical  process  using  mathematical 
equations.  The  implementation  of  those  equations  forms  a 
basis  for  a  simulation  of  the  physical  system.  That  simu¬ 
lation  allows  the  engineer  to  exercise  the  model  of  the 
physical  system  in  a  controlled  manner.  Of  course,  the 
accuracy  of  an  analysis  performed  by  simulation  is  directly 
dependent  on  the  accuracy  of  the  assumed  mathematical  model 
of  the  physical  system. 


1 


2 


Until  recently,  simulation  was  performed  primarily  on 
analog  computers.  Analog  computers  consist  of  a  number  of 
specialized  electronic  devices  that  respond  to  electrical 
Inputs  in  a  manner  analogous  to  the  response  of  physical 
devices  to  physical  inputs.  The  interconnection  of  the 
analog  computer  devices  in  a  way  that  represents  the 
physical  process  leads  to  an  analog  simulation  of  the 
physical  system.  Once  completed  there  will  be  a  one-to-one 
correspondence  between  the  variables  describing  the  physi¬ 
cal  process  and  the  variables  of  the  analog  simulation. 

This  direct  correspondence  makes  the  analog  simulation  an 
excellent  tool  for  studying  the  physical  system  without 
physically  realizing  that  system.  However,  the  initial 
cost,  the  cost  of  maintenance,  the  cost  of  operation  and  an 
inherent  difficulty  in  reprogramming  make  analog  computers 
unattractive. 

As  an  alternative,  modern  digital  computers  offer  easy 
reprogramming  at  relatively  low  costs.  Additionally,  digi¬ 
tal  computers  effectively  solve  multiple  problems  simulta¬ 
neously,  do  not  have  to  be  dedicated  to  simulation  and  pro¬ 
vide  the  ability  to  store  vast  amounts  of  data.  The  com¬ 
bination  of  computational  speed,  programming  flexibility, 
and  sophisticated  numerical  algorithms  for  the  solution  of 
complex  mathematical  equations  has  made  simulation  via 


digital  computers  a  common  practice  over  the  last  fifteen 


years, 


a  Jti 


Many  physical  processes  occur  in  the  environment  that 


can  not  be  precisely  modeled  in  a  deterministic  manner. 
Examples  of  such  processes  are  wind  gusts,  electronic  sen¬ 
sor  noise,  future  target  maneuvers,  the  weather  and  the 
economy.  Although  such  processes  can  not  be  modeled  deter¬ 
ministically,  they  can  be  modeled  as  continuous  random 
processes  based  on  statistical  data.  To  model  a  physical 
process  as  a  random  process,  one  only  needs  an  adequate 
statistical  representation  of  the  physical  process.  This 
representation  can  be  derived  from  previous  observation 
data  or  simply  defined  based  on  knowledge  of  the  physical 
limits  of  the  process.  Although  such  a  derived  or  defined 
model  can  not  precisely  predict  the  behavior  of  the  physi¬ 
cal  system,  it  can  predict  the  statistical  behavior  of  the 
system  or,  in  other  words,  how  the  system  will  behave  "on 
the  average".  Modeling  the  non-deterministic  system  in 
this  manner  will  provide  useful  information,  much  more  so 
than  simply  ignoring  the  noisy  process.  Including  such 
random  models  in  a  simulation  will  certainly  enhance  its 
accuracy  by  making  it  more  faithful  to  the  physical  system. 

fiaadfla  Y.ailablcg 

Before  proceeding  further,  it  will  be  useful  to  define 
and  show  a  number  of  useful  mathematical  properties  related 
to  random  variables.  To  do  this,  let  j.  be  an  n-dimensional 
vector  of  random  variables.  All  of  the  information  known 
about  &  will  be  embodied  in  its  probability  density  func- 


tion  (pdf)  or,  if  the  pdf  is  not  available,  in  its  proba¬ 
bility  distribution  functional)  For  convenience,  assume 


the  pdf  of  &  is  available  and  let  it  be  denoted  by  f^(  • ) . 

f ^(  ■ )  can  be  used  to  compute  the  "expected  value"  of 
some  function  of  &,  where  the  "expected  value"  is  a  numeri 
cal  value  obtained  by  averaging  the  outcomes  of  an  experi¬ 
ment  over  an  ensemble  of  trials.  Specifically,  if  the 
m-dimensional  vector  function  2.(  • )  is  a  function  of  &  such 
that 


2.(0  «  0Cx.(OD 

where  0(0  is  continuous,  then  the  expectation  of  &, 
denoted  as  ElxJ*  is 


E[  2.)  * 


a(i>f*(i>di 


Since  by  definition,  the  expectation  is  an  integral, 
it  is  a  linear  operation.  Therefore,  if 

H.(  • )  ■  A2.(  • )  *  A 0Cx.(  • )  ] 
where  A  is  a  constant  matrix,  then 


E(&)  =  E(AjJ 


fAKPfj/Pdi  «  A  Bipfj 

— •  *  — * 


(Pd£  -  AEI2.1  (1 


Additionally,  if 

!L<  *  >  *  2.1  ( * )  +  i2  ( * )  ■  0j  I  Jt(  • )  1  ♦  J2 1  JL(  * )  1 

then 


El*.]  -  El 2.1 +2.2 J  ■  p®l<P  +  ^(P^jiPdl 


fsitPf^Pdl  +  [^(Pf^Pdl 
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-  ElliJ  +  E l£2J  <2) 

Particular  functions  of  &  are  used  to  characterize 
f^(  • ) .  The  expected  value  of  these  particular  functions 
are  called  the  moments,  or  statistics,  of  x..  Generally,  an 
infinite  number  of  moments  are  needed  to  completely  charac¬ 
terize  f • ) .  Of  special  interest  in  stochastic  analyses 
are  the  first  two  moments  of  &•  The  first  moment  is  called 
the  mean  of  x.  and  the  second  moment  is  called  the  autocorre 
lation  matrix  of 

To  generate  the  first  moment  of  or  the  mean  of  £, 
let  8(&)  ■  The  mean  of  x.»  denoted  as  a,  will  be  compu¬ 


ted  by 


OL  -  E[jl1 


*  f 


<I>d£ 


Note  that  m.  is  not  a  random  variable,  but  rather  a  deter¬ 
ministic  quantity.  This  fact  is  true  for  all  the  statis¬ 
tics  of  a  random  variable. (1] 

To  generate  the  autocorrelation  matrix  of  &,  let 
8(X,)  *  xxT,  where  the  superscr ipt-T  denotes  a  vector  trans¬ 
pose  operation.  The  autocorrelation  matrix  of  &,  denoted 


as  *,  will  be  computed  by 


*  ■  Ef xxTI 


<I)dI 


A  statistic  closely  related  to  the  autocorrelation 
matrix  is  the  covariance  matrix,  denoted  as  P.  Like  the 
autocorrelation  matrix,  the  covariance  matrix  is  a  second 


moment.  P  is  the  second  central  moment  of  &.  The  covari¬ 
ance  matrix  is  defined  as 


P  A  El<*-a>(aL-tt)T]  ■  J  (I  "  tt><£  " 

Using  the  facts  that  El • 1  is  a  linear  operation  and  a 
is  a  deterministic  quantity,  the  relationship  between  *  and 
P  can  be  derived  in  a  straight  forward  manner.  Consider, 

P  -  El  (JL~a)  (X.~a)T3  ■  E[  xxT  -  xmT  -  roxT  +  romT] 

From  equation  (2)  it  is  seen  that 

E  [  xxT -xmT -mxT-f-mmT  ]  ■  Ef  xxT  ]  -  Efxm^]  -  Etn>xT  ]  +  EfmmT  1 
Since  a  is  a  deterministic  vector,  from  equation  (1)  it  is 
seen  that 

e [ smT i  *  Etaia1  *  maT; 

Ef mxT]  *  aE( iT 3  *  mmT: 
and  ElmaT]  ■  mmT. 

Using  these  relationships,  then 

p  -  e[jultj  -  maT  -  *  -  E[*iE[aT]  O) 

Note  that  if  a  *  ®  then  P  ■  *  . 

A  useful  relation  is  statistical  independence.  Two 
random  variables  are  defined  to  be  statistically  indepen¬ 
dent  if  their  joint  probability  density  function  is  equal 
to  the  product  of  their  individual  density  functions. 


Given  two  random  variables,  a  and  v,  then  u  and  v  are 


statistically  independent  if  [2] 

£uv(u,v)  «  fu(u)fv(v) 

where 

fuv(u,v)  4  joint  pdf  of  a  and  v; 
£u(u)  A  pdf  of  a; 
and  fv(v)  =  pdf  of  v. 


Uncorrelated  and  Orthogonal  Random  Variables 


Two  random  variables,  u  and  v,  are  uncorrelated  if  [2] 
E(uv]  ■  E(u)E[v]. 

They  are  orthogonal  if  (2] 

Etuv]  *  0. 


Theorem  1 

If  two  random  variables,  u  and  v,  are  independent, 
they  are  also  uncorrelated. [ 1 ] 

Proof: 


E(uv) 


(u, v)dudv 


(u)du 


( v)dv 


*  E(u]E[v] 


Theorem  2 

If  two  n-dimensional  random  variable  vectors,  &  and  £, 
are  uncorrelated,  then  the  covariance  of  their  sum  equals 


the  sum  of  their  covariances . C 1 1 
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Proof: 


E[((x.-Bx)  +  <  JL-ffly )  H  ( fO*  >  ♦  (JL-ffiy))Tl  “ 

Et  (&-Bx>  (X.-&x>Tl  +  E [  ( l-m* >  ( X.-ffiy ) T 1  +  Et  ( 3£_— ffly )  (X.-Bat)TJ 


E[  (2£.-ay)  (Y.-ffly)Tl 


Enx-B*MX-ffl*>T3  *  Px; 

Et  (JL-By)  (X.-m,y)Tl  “  Pyi 

Et  (s.-blx)  (3L-ffly)Tl  -  EtiY.TJ  "  OatEt^J  -  EtjtJOy  +  m.x&y  “ 


Et  ( 3L“lBy H i-m* ) T 1  *  Et3UtTl  “  myEtiT]  -  EtYJaJ  *  %Ex  =  0* 


Therefore, 


EKtlrl#*)  ♦  ( JL-ffly )  >  ( ( X.-B* >  +  (2L“ly)>TJ  *  Px  +  Py 


Theorem  2  can  easily  be  extended  to  provide  the  relation¬ 


ship  that  if  the  random  vector. 


x  -  2$  u , 
i-11 


then 


PX(D 


Likewise,  from  equation  (1)  and  Theorem  2,  if 


then 


E  “  f?iAi*i 
Pxd)  “  ^jj^iPui*! 


T  141 


A  process  that  contains  an  element  of  chance  is  for¬ 


mally  called  a  stochastic  process.  A  stochastic  process 


will  be  described  by  some  function  that  contains  one  or 


more  random  variables.  The  statistics  of  the  stochastic 


process  trill  be  a  function  of  not  only  the  statistics  of 
its  random  variables  but  also  the  describing  function  it¬ 
self.  Some  examples  of  stochastic  processes  are 

X(t)  ■  A  sin(wt)  where  u  is  a  random  variable  with  a 

known  pdf. 

?(t)  =  some  noiselike  signal  with  no  deterministic 
structure . 

x(t)  »  w(t)  where  w(t)  is  a  random  varible  with  a 

uniform  distribution. 

Ingredients  of  Stochastic  Simulation 
Analog  Simulation 

In  order  to  include  a  random  process  in  an  analog  sim 
ulation,  one  needs,  in  addition  to  the  mathematical  model 
of  the  process,  an  electronic  device  that  will  generate 
wide-bandwidth  noise.  Hide-bandwidth  noise  is  noise  that 
contains  power  components  over  a  very  large  range  of  fre¬ 
quencies.  Ideally,  the  noise  generator  should  produce 
noise  that  has  a  frequency  spectrum  of  constant  amplitude 
over  an  infinite  range  of  frequencies.  Such  noise  is  said 
to  have  an  infinite  bandwidth  with  a  constant  spectral 
density  function  and  is  termed  "white  noise". 

The  desire  to  produce  white  noise  comes  from  the  fact 
that  the  mathematics  involved  with  analyzing  systems  con¬ 
taining  white  noise  is  greatly  simplified.  For  instance, 
the  autocorrelation  function,  *(T),  for  a  white  noise 


vector  is  of  the  form  (3] 


where 


*<T)  «  Q  «tt) 


f(?)  s  Dirac  delta  function 
and  Q  is  the  spectral  density  matrix  of  the 
noise 

The  Dirac  delta  function  is  a  defined  mathematical  function 
that  has  infinite  amplitude,  zero  width,  and  occurs  at 


T  *  0. 


As  seen  from  its  autocorrelation  function,  white  noise 
is  not  correlated  in  time.  In  other  words,  knowing  the 
value  of  the  noise  at  any  time  provides  no  information 
concerning  the  value  of  the  noise  at  any  other  time.  This 
uncorrelatedness-in-time  property  as  well  as  the  constant 
spectral  amplitude  property,  vastly  simplifies  the  mathema¬ 
tics  involved  with  processing  the  noise  statistics.  Thus, 
from  a  mathematical  sense,  a  white  noise  model  is  a  very 
attractive  model.  However,  white  noise  is  physically  not 
possible.  If  this  isn't  obvious  from  the  fact  that  it 
contains  power  components  at  an  infinite  number  of  frequen¬ 
cies,  consider  the  argument  that  white  noise  changes  it's 
value  by  an  infinite  amount  in  zero  time.  One  might  ask 
how  such  a  noise  model  could  be  justifiably  used  in  a 
simulation  of  a  real  system.  There  are  two  justifications. 

First,  if  a  physical  noise  generator  has  a  fairly  flav 
spectral  density  over  a  range  of  frequencies  that  is  much 
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greater  than  the  bandwidth  of  the  system  that  Is  driven  by 
the  noise  source,  then  the  effect  of  the  bandlimited  noise 
on  the  system  is  approximately  the  same  as  if  it  were  driv¬ 
en  by  white  noise.  Figure  1  demonstrates  this  point  graph¬ 
ically.  [1]  Second,  virtually  any  spectral  density  function 
can  be  shaped  from  a  white  noise  spectral  density  function 
by  processing  the  white  noise  through  a  shaping  filter. [1] 


Therefore,  the  model  of  the  true  noise  source  can  be  simu¬ 
lated  by  a  white  noise  source  driving  a  shaping  filter. 

The  white  noise  source  would  in  reality  be  a  wide-bandwidth 
noise  generator  that  has  a  bandwidth  much  greater  than  the 
bandwidth  of  the  shaping  filter.  Since  an  analog  simula¬ 
tion  uses  continuous  analog  devices  and  models  to  simulate 
the  continuous  noise  process,  there  is  a  one-to-one  corre¬ 
spondence  between  the  statistics  of  the  simulated  noise 
process  and  the  statistics  of  the  actual  noise  process. 


This  one-to-one  correspondence  is  indeed  a  nice  relation¬ 
ship  because  it  precisely  defines  the  required  noise  sta¬ 
tistics  to  be  used  in  the  simulation  in  order  to  accurately 
model  the  assumed  physical  noise  process.  Obviously,  it 
would  be  equally  attractive  if  this  one-to-one  correspon¬ 
dence  held  for  a  digital  simulation.  Unfortunately,  it 
does  not. 

Digital  Simulation 

To  simulate  a  continuous  random  process  in  a  digital 
simulation,  one  needs  the  same  functional  ingredients  as  is 
required  in  an  analog  simulation;  that  is,  a  noise  genera¬ 
tor,  a  shaping  filter,  and  knowledge  of  the  relationship 
between  the  statistics  of  the  continuous  noise  process  and 
the  statistics  of  the  digital  simulation  model.  Unlike  an 
analog  computer,  a  digital  computer  is  a  device  that  is 
only  capable  of  producing  a  finite  number  of  conditions  or 
states.  Because  a  digital  computer  is  a  finite  state 
machine,  it  is  not  possible  to  generate  a  truly  continuous 
random  process.  The  best  that  can  be  achieved  is  the 
generation  of  a  repeatable  sequence  of  numbers  that  within 
the  finite  window  of  a  sample  space  appears  to  be  random. 
However,  if  the  sample  is  large  enough,  this  "random" 
sequence  of  numbers  is  usually  adequate  for  simulating 


continuous  noise. 
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Methods  for  generating  random  sequences  on  a  digital 
computer  (commonly  called  random  number  generators)  are 
bountiful.  Many  of  these  random  number  generators  provide 
excellent  and  predictable  statistical  properties.  The  most 
commonly  used  random  number  generators  produce  a  random 
sequence  that  has  a  uniform  distribution  over  some  range  of 
numbers.  This  uniform  distribution  can  be  shaped  into  most 
desired  noise  distributions  through  the  use  of  shaping 
filter  algorithms  in  much  the  same  way  as  is  done  in  analog 
computers . 

The  final  ingredient  needed  to  simulate  continuous 
noise  in  a  digital  simulation  is  the  knowledge  of  the  rela¬ 
tionship  between  the  statistics  of  continuous  noise  process 
and  the  digital  or  discrete  noise  process.  Unlike  in  ana¬ 
log  simulation,  there  is  not  a  one-to-one  correspondence 
between  these  statistical  properties.  In  fact,  this  dis¬ 
sertation  will  show  that  the  relationship  is  complex  and  a 
factor  of  not  only  the  continuous  noise  statistics  and  the 
dynamics  of  the  system,  but  also  of  the  numerical  algor¬ 
ithms  used  to  solve  the  equations  of  motion. 

The  simplest  situation  in  which  a  relationship  would 
be  required  is  the  case  when  one  is  attempting  to  simulate 
a  continuously  available  measurement  of  a  state  that  is 
corrupted  by  additive  white  noise.  This  situation  is 
described  mathematically  by 

ze(t)  -  x ( t )  +  ve(t) 
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where 

zc(t)  *  continuous  measurement  as  a 
function  of  time,  t; 

x(t)  *  a  deterministic  function  or  state; 
and  vc(t)  *  a  continuous  white  noise  function. 
The  discrete  measurement  equation  for  this  case  is 

Zd(ti)  -  x(tA)  +  Vjjtt^) 

where 

Zd<ti)  =  discrete  measurement  taken  at 
time  equal  to  tj; 

x(ti)  =  the  value  of  x(t)  at  t  -  t^; 
and  v^t^)  “  a  discrete  noise  term  applied  at 

t  *  t*. 


In  a  simulation  that  contains  this  type  of  model,  it 
is  necessary  that  the  statistics  of  z^lt^)  equal  the 
statistics  of  zc(t)  evaluated  at  t  *  t^,  for  all  t^.  Since 
x(t)  is  deterministic  and  the  noise  is  additive,  this  re¬ 
quirement  will  be  met  if  the  statistics  of  v^Ct^)  equal  the 
statistics  of  vc(t)  evaluated  at  t  *  t*,  for  all  tA. 

Maybeck  (1)  develops  the  relationship  to  meet  this  require¬ 
ment.  A  summary  of  Maybeck’s  development  is  given  below. 

Assume  the  vc(t)  is  a  zero-mean  Gaussian  noise  with 

Etvc(t)vc(t+T) ]  =  RC«(T) 

where 

E(>]  denotes  the  statistical  expectation; 

Rc  *  spectral  density  of  vc(t); 


and 


« (T)  denotes  the  Dirac  delta  function. 


If  v^Ct^)  is  a  zero-mean  Gaussian  sequence  with 

EIvd2(ti)l  -  R^t^  *  Rc/Ati  (6) 

where  AtA  A  t^+i  -  t^ 

and  E( vd( t^) vd( t j ) )  «  0  for  i  *  j 

then,  in  the  limit  as  At^  -*  0,  equation  (6)  will  reach  the 
strength  of  R(*<(0).  Hence,  the  desired  relationship 
between  the  covariance  of  the  random  number  generator  and 
the  spectral  density  of  the  continuous  noise  process  is 
given  by  equation  (6).  Note  that  even  in  this  simple  case 
that  includes  no  dynamics,  there  is  is  not  a  one-to-one 
statistical  correspondence.  Discretization  has  introduced 
a  functional  relationship  between  the  statistics  of  the 
simulated  noise  process  and  the  noise  process  of  the  simu¬ 
lator.  The  relationship  given  by  equation  (6)  is  the  sim¬ 
plest  relationship  that  will  be  developed  in  this  study. 
When  the  process  equations  contain  dynamics,  the  functional 
relationships  will  become  more  complex  and,  without  simpli¬ 
fying  assumptions,  will  also  be  explicit  functions  of  the 
dynamics . 

Dissertation  Brief 

Specifically,  the  class  of  problems  that  will  be  ad¬ 
dressed  in  this  dissertation  are  linear  first-order  differ¬ 
ential  systems  that  are  driven  by  Gaussian  white  noise. 
Relationships  will  be  derived  for  the  accurate  solution  of 
this  class  of  problems  in  digital  simulations  that  use  a 


number  o£  popular  integration  methods.  It  will  be  shown 
that  each  integration  method  yields  a  different  relation¬ 
ship  for  the  noise  statistics.  It  will  further  be  shown 
that  a  number  of  popular  integration  methods  that  are  often 
preferred  in  deterministic  simulations  are  impractical  for 
simulations  of  stochastic  systems. 

The  general  problem  will  be  described  and  analytically 
developed  in  Chapter  II.  The  necessary  functional  relation 
ships  between  the  continuous  noise  process  statistics  and 
the  discrete  noise  statistics  for  the  integrators  of 
interest  will  be  derived  in  Chapter  III.  An  error  analysis 
for  the  results  obtained  in  Chapter  III  will  be  developed 
in  Chapter  IV.  The  results  from  this  analysis  will  define 
the  limits  of  validity  for  a  simplifying  assumption  made  in 
Chapter  II.  The  analytical  findings  of  this  research  will 
be  demonstrated  through  simulation,  comparing  the  numerical 
results  to  known  analytical  solutions.  The  simulation 
results  will  be  given  in  Chapter  V.  The  research  findings 
will  be  summarized  and  recommendations  will  be  made  in 
Chapter  VI. 


II.  PROBLEM  DESCRIPTION 

This  research  is  applied  to  the  analysis  of  problems 
belonging  to  the  class  of  linear  time-invariant  stochastic 
differential  systems.  Mathematically,  this  class  of  prob¬ 
lems  is  described  by 

i(t)  -  A*(t)  +  Bua.(t)  +  Bw»L(t)  (7) 

where 

X.(t)  =  n-dimensional  state  vector; 


!i<t)  =  p-dimensional  deterministic  input  vector; 
St( t )  =  m-dimensional  random  input  vector; 

A  =  nxn  system  matrix  of  constant  elements; 

Bu  =  nxp  input  matrix  of  constant  elements; 

and  Bw  =  nxm  input  matrix  of  constant  elements. 

For  problems  in  the  form  of  equation  (7),  the  state  vector 
can  be  described  by  the  summation  of  two  n-dimensional  vec¬ 
tors,  one  purely  deterministic  and  the  other  purely 
stochastic.  Thus, 

S.(t)=j^,(t)+xw(t)  (8) 

where 

^(t)  =  deterministic  state  vector 
and  X«(t)  =  stochastic  state  vector. 
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Substituting  equation  (8)  into  equation  (7)  yields 


iu(t)  ♦  x*(t)  *  CAiu(t)  +  Bua.(t)3  +  EAj*(t)  +  Bwa<t)D  (9) 


Clearly  from  equation  (9)  the  problem  of  analyzing  the  sys¬ 


tem  can  be  separated  into  a  deterministic  analysis  problem 


and  a  stochastic  analysis  problem.  Since  the  stochastic 


analysis  problem  is  the  one  of  interest,  it  will  be 


assumed,  without  loss  of  generality,  that  &(t)  =0  and 


^(0)  =  0.  By  making  these  simplifying  assumptions  and 


dropping  the  subscripts,  one  finds  the  stochastic  system 


to  be  described  by 


i<t)  *  Aj.(t)  +  Bg.(t) 


(10) 


Since  &(t)  is  a  non-deterministic  function,  equation 


(10)  can  not  be  solved  explicitly.  The  best  that  one  can 


do  is  solve  for  the  statistics  of  &(t),  which  are  determin¬ 


istic  quantities.  If  jf.(t)  is  Gaussian,  the  first-order  and 


second-order  statistics,  mean  and  covariance,  fully 


describe  the  statistical  properties  of  g.(t).[l,2J  Further, 


since  the  system  is  linear,  the  statistics  of  &(t)  will 


also  be  Gaussian. [1]  Thus,  one  only  needs  to  determine  the 


mean  and  covariance  of  &(t)  to  determine  the  complete  sta¬ 


tistics  of  x.(t).  If  g.(t)  is  non-Gaussian  then  higher  order 


statistics  will  be  needed  for  a  complete  statistical  de¬ 


scription.  ( 2]  However,  even  for  non-Gaussian  distributions. 


the  system  analyst  is  often  only  interested  in  determining 


the  first  two  statistical  orders  due  to  the  physical 


meaning  o£  those  statistics.  For  convenience,  it  will  be 
assumed  that  g.(t)  has  a  Gaussian  distribution. 


Determination  of  the  Statistics  of  the  State  Vector 
Since  equation  (10)  is  a  stochastic  differential  equa¬ 
tion,  normal  solution  techniques  can  not  be  used  to  deter¬ 
mine  &(t).  Another  approach  will  be  used  in  order  to  find 
the  form  of  &(t)  that  will  allow  the  determination  of  the 
statistics  of  &(t).  To  this  end,  define  a  random  variable 
vector,  £.(t),  such  that 

£(t)  =  f  £<*>  dt 

Jg 

where  h.(T)  is  a  zero-mean  Gaussian  white  noise  vector  with 
a  spectral  density  matrix  equal  to  Q.  £.(t)  is  called  a 
Weiner  or  Brownian-motion  process . ( 1, 2, 3 ]  The  statistics 
of  g.(t)  can  easily  be  determined.  The  mean  or  average  of 
Brownian-motion  is 


E(&(  t)  ] 


-  e[  f«x> ar  ]  ■  I; 


Era.m  i  dt  =  0 


The  autocorrelation  matrix  of  §.(t)  is 

E[g.(t)g.T(t)  ]  *  e[  |  3.(1)  dX  |  jLT(«r)  d<r  j 


ff- 


E(£(X)wT(r) ]  dtd<r 
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As  presented  in  Chapter  I,  the  E[H.(T)wT(cr)  1  =  Q*(T-<r)  ,  there¬ 
fore 

Q*<T-0  dTdc  =  I  Q  d<r  =  Qt 
*»0 

Using  this  defined  function,  £(t),  equation  (10)  can  be 
rewritten  in  another  form  that  can  be  solved  for  the  statis¬ 
tics  of  &(t).  Since  equation  (10)  is  a  linear  differential 


E(g.(t)g.T(t)  ] 


-IT 


I 

fl 


equation  in  time,  it  can  be  rewritten  as  [1] 

di(t)  =  A&(  t  )dt  +  Bdg.(  t ) 

This  equation  can  be  integrated  to  yield 

*(t)  =  l(t0)  +  f  Ax.(T)  d*t  +  f  Bdfc(T) 

Jtg  J tg 

The  problem  at  hand  is  to  determine  a  stochastic  process 
l(  • )  that  satisfies  this  integral  equation.  In  section  4.8 
of  Maybeck  [1],  it  is  shown  that  the  following  equation  is 
such  a  process. 

l(t)  =  $( t, tg )  &(  tg )  +  f  5(t,T)Bdg.m  (11) 

Jtg 

where  $(t,tg)  is  the  state  transition  matrix  that  satisfies 
5(t,tg)  *  A5(t,tg)  and  3(tg,tg)  =1,  where  I  is  the  ident¬ 
ity  matrix  and  A  is  the  system  matrix  in  equation  (10). 


Note  that  x.(t)  is  a  stochastic  process,  thus,  its  solution 
is  non-deterministic;  however,  using  the  knowledge  of  the 
statistics  of  £(t),  the  statistics  of  j,(t)  can  be  deter¬ 
mined. 


(12) 


The  mean  of  j,(t),  denoted  E(x.(t)l,  is  described  by 
E(i<t)l  *  5(t,t0)E{jL(t0)  ]  +  E[  |  5<t,T)Bdg.<T)  J 

where  5(t,t0)  =  fA(t-t0) 

-  I  +  A*(t-t0)  +  &2*(t-t0)2  +  .  .  . 

2! 

and  t  and  t0  denote  arbitrary  values  of  time.  In  equation 
(12),  I  is  an  nxn  identity  matrix.  Since  E(g.(t)]  =  0,  the 
stochastic  integral  is  also  zero-mean  and 

Ett(t)]  =  2( t,  t0)E[&( t0) ]  (13) 


The  autocorrelation  matrix  of  j.(t)  can  be  determined 
using  equation  (11)  and  equation  (2)  to  write  Etjt.(  t)£T(  t)  ] 
as  the  sum  of  four  separate  expectations.  Using  the  defi¬ 
nition  of  &(t),  x.(t0)  is  independent  of  &(t).  From  Theorem 
1  and  the  fact  that  the  mean  of  the  stochastic  integral  is 
zero,  the  expected  value  of  two  cross  terms  is  zero.  For 
example, 

ft  at 

5(t,T)Bd6(T)xT(to)l  =  eT  S ( t , T ) BdB ( T ) 1 E ( xT ( tn ) ]  =  0 


r  f4  -i 

,  r  ffc  i 

5(t,T)Bd£(T)i.T(t0) 

=  E  2(  t ,  T ) BdjB(  T ) 

L  J*0  J 

1  L  Jt0  J 

Therefore,  the  autocorrelation  matrix  of  &(t)  will  be 
E[i(t)£T(tH  =  5(t,t0)EIi(t0)s.T(t0)  ]2T(t,t0) 


r«t. 

Jtfl 


T  )BQ(  *t  )BT5T(  t ,  T  )dT 


The  covariance  of  £(t),  denoted  as  P(t),  can  be  de¬ 
rived  directly  from  the  autocorrelation  matrix  of  z.(t)  by 
the  relationship  given  in  equation  (3).  Specifically, 
Etl<t)iT(t)  ]  =  P(t)  +  E [ x.( t )  ] E [ xJ ( t )  ] 
and  EtjL(ta)x.T(tfl)  J  =  P(t0)  +  E[*(  t0)  ]E[jLT(t0)  ] 

Evaluating  the  autocorrelation  matrix  of  &(t)  using  these 
relationships  and  equation  (13),  yields 
P(t)  =  5<t,t0)P(t0)5T(t,t0) 


[\<t. 

Jt0 


T)BQ(T)BT5T(t,T)  dt 


(14) 


Notice  that  equation  (14)  is  a  deterministic  equation  as  is 
equation  (13). 


A  stochastic  analysis  of  this  system  involves  solving 
for  the  statistics  of  &(t);  specifically,  for  the  mean  and 
covariance  of  &(t).  For  low-order  linear  systems,  equation 
(13)  and  equation  (14)  can  be  solved  analytically;  but  for 
large-order  linear  systems  and  of  course  nonlinear  systems, 
simulationr  die  normaiy  used  to  find  the  statistics  of 
&(t).  In  a  digital  simulation,  a  random  number  generator 
is  used  to  generate  a  random  sequence,  that  will  simu¬ 

late  K.(t).  By  repetitively  running  the  simulation  using  a 
large  number  of  sample  sequences  from  the  random  number 
generator  and  statistically  averaging  the  results,  an 
approximation  of  E[j.(t)]  and  P(t)  can  be  generated.  This 
analysis  method  is  commonly  called  a  Monte  Carlo  analysis. 
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In  a  digital  simulation,  g(t)  will  be  represented  by  a 
discrete  state  vector  go  and  the  propagation  o£  g.(t)  will 
be  approximated  by  the  propagation  of  gjj  using  discrete 
difference  equations.  In  order  to  insure  a  faithful  simu¬ 
lation,  the  statistics  of  g^tt^)  must  be  approximately 
equal  to  the  statistics  of  g(t)  evaluated  at  t  *  for  all 
t  over  the  range  of  time  of  interest.  This  can  be  accom¬ 
plished  on  an  incremental  basis  assuming  that  the  statis¬ 
tics  of  g£(tg)  equal  the  statistics  of  g,(tg)  where  t0  is 
some  arbitrary  initial  time.  By  insuring  that  the  propaga¬ 
tion  of  the  statistics  of  g£  approximates  the  propagation 
of  the  statistics  of  g.  from  one  update  time  to  the  next  and 
assuming  that  accumulated  errors  remain  small,  then  the  re¬ 
quirement  for  a  faithful  simulation  will  be  realized. 

Using  equation  (13),  the  equation  that  describes  the 
propagation  of  the  mean  of  g.(t)  from  t  ■  tj  to  t  *  ti  +  l 
where  t^+i  *  t^  +  h  will  be 

ECX.C ti+i )  I  ■  5(h,0)ECg.(ti)  ]  (15) 

Likewise,  from  equation  (14),  the  covariance  of  g.(t)  will 
be  propagated  from  t  *  ti  to  t  =  ti+1  by 
P(ti+1)  =  5(h,0)P(ti)5T(h,0) 

fh 

+  5(h,T)BQBT5T(h,T)  dT  (16) 

J  0 


In  the  digital  simulation,  the  discrete  state  £q  would 
be  propagated  from  t ±  to  t^  +  i  *  (t^  +  h)  by  [1] 


XO(ti  +  l)  =  S^h^)*^  tA)  +  BjjSiotti) 


(17 


where 

5d(»,*)  is  a  discrete  state  transition  matrix 
and  Bq  is  a  discrete  input  matrix. 

5DC'->  will  not  only  be  a  function  of  the  A  matrix  in 
equation  (10)  and  Bq  will  not  only  be  a  function  of  the  A 
and  B  matrices  in  equation  (10),  but  these  two  functions 
will  also  be  dependent  on  the  integration  method  used  in 
the  simulation. 

Determination  of  the  Discrete  State  Statistics 
Mean  of  the  Discrete  State  Vector 

If  the  mean  of  stg  is  defined  to  be  zero  for  all  ti# 
then,  from  equation  (17),  it  is  seen  that  the  mean  of  xg 
will  be  propagated  from  t^  to  tj+i  by 

EtxoCti+i)]  =  5D(h,0)Etso(ti) ]  (18) 

Second  HaasuL  aJL  geests,  Stats  Jtegi&L 

The  autocorrelation  matrix  of  xq  can  be  derived  direct¬ 
ly  from  equation  (17)  by  determining  Etxj)( t^+^)  ] . 
Since  the  discrete  state  xo^i*  can  not  be  influenced  by  the 
input  Xo(tj.)  (present  state  conditions  can  only  be  influ¬ 
enced  by  past  inputs  and  initial  state  conditions),  Xo^i* 
is  independent  of  j^ft^).  Therefore,  by  Theorem  1  and  the 
fact  that  Elg^t^)]  =  0,  the  expected  value  of  the  two 
cross  terms  will  be  zero.  Thus, 

EIXo(ti)*S(  ti)  3  =  E[jj0(ti)xJ(ti)  ]  =  0 


EtioCti  +  iJjJtti+i)  ]  =  SQOuajEtXDCtiJxgUi)  )S$(h,0) 

+  bdqdb5 

where  0q  ■  E Cj£i)(  t^ ) w£(  t^ )  3  for  all  t^. 

The  discrete  covariance  matrix,  Pq(>),  can  be  deter¬ 
mined  from  the  autocorrelation  matrix  of  by  using  the 
following  relationships 

Etlo(ti+l)x!(ti+l) 1  "  pD(ti+l>  +  EtXo(ti+1) ]EtiT(ti+1) ] 

E[x0(ti)xJ(ti)  ]  =  PD(ti)  +  EIxoUi)  JEt^tti)  ] 
and  E[£D<ti+1)]  =  5D(h,0)Etxo( tL) }  =  E[ 5D( h, 0 )xD( tA ) ] . 

Using  these  relationships, 

PD(ti+1)  =  5D(h'0)pD<ti)5|5(h'0)  +  bdQdb5  (19) 


Development  of  Conditions  for  a  Faithful  Simulation 
By  direct  comparison  of  equation  (18)  to  equation 
(15),  it  is  clear  that  in  order  for  the  propagation  of  the 
mean  of  the  discrete  state  vector  to  approximate  the  propa¬ 
gation  of  the  mean  of  the  continuous  state  vector,  5D(h,0) 
must  be  approximately  equal  to  5(h,0).  This  is  the  same 
requirement  needed  to  insure  digital  simulations  of  deter¬ 
ministic  systems  are  faithful.  Thus,  it  will  be  assumed 
that  the  numerical  integrator  used  in  the  digital  simula¬ 
tion  will  insure  that 

3D(h,0)  =  2(h,0) 

Making  this  assumption,  a  comparison  of  equation  (19)  to 
equation  (16)  yields  that  the  propagation  of  the  discrete 
covariance  will  be  approximately  equal  to  the  propagation 


of  the  continuous  covariance  if 

fh 

BDQdb8  =  5<h,T)BQBT2T(h,T)  dT  (20) 

"  0 

Selection  of  an  integration  method  that  will  satisfy 
3[)(h,0)  =  5(h,0)  will  insure  that  EEs^Ct^)]  approximates 
E[S.(t)]  evaluated  at  t  =  t^.  However,  additional  require¬ 
ments  are  needed  to  ensure  that  the  second-order  statistics 
are  faithfully  simulated.  Thus,  equation  (20)  describes 
the  essence  of  the  problem.  If  3(»,*)  is  known  then  the 
right  side  of  equation  (20)  can  be  solved  either  in  closed 
form  or  numerically.  However,  in  general  £(•,•)  will  not 
be  known.  Of  course,  the  general  case  is  the  one  of  inter¬ 
est.  What  is  known,  in  general,  is  that  as  h  approaches 
zero,  3(h,0)  approaches  the  identity  matrix.  In  previous 
works  of  this  type,  this  fact  was  used  routinely  to  approx¬ 
imate  the  state  transition  matrix  with  the  identity  matrix 
in  order  to  evaluate  the  right  side  of  equation  (20). 

[1,3,51  At  the  outset,  the  analyses  developed  in  this  work 
will  also  use  this  relationship.  However,  because  of  the 
practical  significance  of  this  convenient  assumption,  later 
in  the  error  analyses  section.  Chapter  IV,  this  assumption 
will  be  fully  investigated. 

If  h  is  selected  reasonably  small  such  that 

h*|\nax|  «  1 

where  |^max|  is  the  magnitude  of  the  largest  eigenvalue  of 


A,  then 
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£(h,0)  s  I 


(21) 


Using  equation  (21),  evaluation  o£  equation  (20)  yields 


BqQqBJ  =  hBQBT 


(22) 


In  a  digital  simulation  used  to  perform  a  Monte  Carlo 
analysis  of  this  type  system,  it  is  necessary  to  select  an 
h  that  will  satisfy  equation  (21)  and  it  is  necessary  to 
adjust  the  covariance  of  the  random  number  generator  such 
that  equation  (22)  is  satified.  It  will  be  shown  in  Chap¬ 
ter  IV  of  this  dissertation  that  the  need  for  an  exception¬ 
ally  small  h  imposed  by  equation  (21)  can  be  relaxed  for 
certain  integration  methods. 

One  additional  problem  is  that  Bq,  the  discrete  input 
matrix,  is  not  known.  Bq  will  be  a  function  of  the  inte¬ 
gration  method  used  in  the  simulation.  For  each  integrator 
of  interest,  Bq  will  be  determined  in  order  to  derive  the 
necessary  conditions  that  will  satisfy  equation  (22). 
Griffith  (4)  showed  that  if  an  Euler  integrator,  which  will 
be  defined  in  Chapter  III,  is  used  in  the  simulation,  then 


Bq  3  hB 


(23) 


He  then  concluded  using  equations  (22)  and  (23)  that  for  a 


simulation  employing  an  Euler  integrator 

QD  =  Q/h 


(24) 


Griffith  then  applied  his  analytical  findings  to  a  problem 
that  used  an  Adams -Bashf or th  integrator  (also  to  be  defined 
in  Chapter  III).  Griffith's  numerical  results  seemed  to 
show  that  equation  (24)  was  also  a  valid  relationship  for 
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the  Adams-Bashforth  method.  It  will  be  shown  analytically 
in  this  work  that  Griffith's  numerical  findings  are  valid 
only  under  very  restricted  conditions.  Further,  it  will  be 
shown  that  each  for  integration  method  there  exists  a 
unique  function  relating  the  continuous  input  noise  statis¬ 
tics  to  the  random  number  generator's  noise  statistics. 
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III.  INTEGRATOR  ANALYSES 


The  numerical  algorithms  that  will  be  considered  in 
this  study  are  designed  to  approximate  the  solution  to  or¬ 
dinary  differential  equations  and  can  be  divided  into  two 
major  categories;  fixed-step  methods  and  multistep  methods. 
Fixed-step  methods  utilize  derivative  information  that  is 
computed  or  approximated  within  the  current  integration 
step  interval,  whereas  multistep  methods  utilize  derivative 
information  accumulated  over  a  number  of  steps.  There  are 
advantages  and  disadvantages  to  each  type  of  integration 
method.  The  type  of  algorithm  that  an  analyst  chooses  to 
use  in  a  given  situation  is  usually  the  one  that  he  feels 
will  provide  accurate  and  timely  solutions.  Many  integra¬ 
tion  techniques  have  been  studied  in  detail  to  assist  the 
analyst  in  making  an  intelligent  choice.  However,  most  of 
these  studies  have  been  performed  with  application  to  de¬ 
terministic  systems  and  not  to  stochastic  systems.  This 
dissertation  shows  that  when  applied  to  solving  stochastic 
differential  equations,  each  integration  method  will  affect 
the  relationship  of  the  noise  generator  statistics  to  the 
continuous  noise  process  statistics.  In  this  chapter,  that 
relationship  between  the  statistics  will  be  developed  for  a 
number  of  widely  used  integration  algorithms. 


Given  a  differential  system  of  the  form 


i  *  £.(  &,  t )  (25) 

a  general  Runge-Kutta  method  approximates  the  solution 
at  time  t  *  (i+l)h,  denoted  as  x^+1,  by  (5] 

S-i+1  *  *i  +  (alli.l  +  «2&.2  +  •••  +  «n!in)  (26) 

where 

Jt!  *  b*£(Si  ,  ih)  (27) 

and 

tj  *  h*t(j.^  +  aj!tj_i,ih  +  bjh)  for  j  =  2,n  (28) 

where  ocj,  <Xj ,  aj,  bj  are  dependent  on  the  particular  Runge- 
Kutta  method  and  £(•#•)  is  defined  in  equation  (25). 

Evilcr  Inteqratai 

An  Euler  integrator  can  be  classified  as  a  Runge-Kutta 
integrator  of  order  one.  For  the  system  of  interest, 

f(JLi*ih)  =  Aj.i  +  Bjti  (29) 

Applying  equation  (29)  to  equations  (26)  -  (28)  with  n  -  1, 
the  Euler  integrator  will  approximate  the  solution  of  j.(t) 
by 

A.i  +  1  “  (I  +  aihA)^  +  c^hBjUn 
where  wqj  is  a  discrete  noise  input  vector  at  the  evalua¬ 
tion  of  f(*,»).  Note  that  for  =  1, 

(I  +  «xhA)  =  (I  +  hA)  =  5D(h,0) 
where  3{}(h,0)  is  a  first-order  truncation  of  3(h,0)  given 
by  equation  (12).  Therefore,  the  update  equation  for  the 


Euler  Integrator  can  be  written  as 

JLi+1  “  , 0 ) Jti  +  aihBjfoi  (30) 

The  discrete  autocorrelation  matrix  for  an  Euler 
integrator  can  be  determined  directly  from  equation  (30). 

Et*i+l*i+lJ  s  5D(h,0)E[x.ixIl5g(h,0)  +  a1hBE(!Dixf]lg(h,0) 
+  5D(h,0)E[iiS£5il°‘ihBT  +  (a1h)2BE[j£oiHjl)BT 
Using  the  relationship  that  the  statistics  of  the  states  at 
present  time  are  independent  of  the  statistics  of  the  pres¬ 
ent  inputs,  and  the  fact  that  E C StDi  1  =  0  for  all  i,  results 
in 

EtKDl*l]  =  EtHoilEt  if]  =  0 
and  EIx^hJiI  »  ECaLi hSi )  =  0 

Therefore, 

E[JLi+1xJ+1]  =  5D(h,0)E(iiif  15jth,0)  +  (a1h)2BE(HoiHJ51)BT 
From  the  autocorrelation  matrix,  the  covariance  matrix  can 
be  determined  using  the  relationships 

EIJLi+nI+lJ  ■  PD(ti  +  l>  "  E[ii  +  1)E[i|+1) 

*  pD(ti  +  l}  "  5D(h'0>Et)Li]E[iJ]5j(h,0) 
and  =  PD(ti)  -  EtijJEtjJ] 

Doing  so  results  in 

PD(ti+i)  =  5D(h,0)PD(ti)5j(h,0)  +  «fh2BQDBT  (31) 

where  Qq  *  E  [  StOlStSi  3  -  Recognizing  that  equation  (31)  is  in 
the  form  of  equation  (19),  provides  the  relationship 

BDQDB5  *  afh2BQDBT  (32) 

Equating  equation  (32)  to  equation  (22)  yields  the  statis¬ 
tical  relationship 
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Qq  *  Q/(aJh) 


(33) 


For  =  1,  equation  (33)  is  the  same  result  that  Griffith 
derived.  (41  Note  that  equation  (33)  is  independent  of  the 
system  dynamics.  This  is  a  nice  relationship  because  it 
means  that  equation  (33)  can  be  applied  to  any  digital  simu¬ 
lation  that  uses  an  Euler  integrator. 


A  second-order  Runge-Kutta  integrator  approximates  the 
solution  to  &(t)  at  t  *  (i+l)h  by  the  folowing  algorithm. 


where 


Z.i  +  1  *  li  +  <«l!i.l  +  «2&-2> 
fc.1  -  h£(ii,ih)  »  h(Aj.i  +  Bsjoi) 

1 *  h^^+ailtl*  (  i+bi)h)  ■  h(  Ali^a]^)  +  Bg^) 
*  [hA  +  ajUiA)2!^  +  alh2ft®H01  +  hBttD2 


(34) 


(35) 


(36) 


where  gQj  and  kq2  ar®  the  discrete  noise  inputs  from  the 
random  number  generator  at  the  respective  evaluations  of 
£(*,*).  Evaluating  equation  (34)  using  (35)  and  (36)  and 
combining  terms  yields 


X.i+1  “  Cl  +  («i+«2)hA  +  a1<X2(hA)2]jti 
+  h(<x^I  +  a^oc2hA )  Bgoi  +  «2hB!£02 


(37) 


The  most  commonly  used  2nd  order  Runge-Kutta  method  defines 
«1  ■  «2  *  ®*5  and  a^  ®  bj_  »  1.  With  these  defined  constants 
note  that 

Cl  +  ( aj-f«2 )hA  +  aj^thA)2]  *  Cl  ♦  hA  +  0.5(hA)2]  =  5D(h,0) 
where  3p(h,0)  is  a  second-order  truncation  of  3(h,0)  given 
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by  equation  (12).  Therefore,  equation  (37)  can  be  rewrit¬ 
ten  as 

*d+l  “  5[)(h,0)x.i  +  BDlH01  +  bD2S02  <38) 

where 

®D1  —  h(  a^I  +  a^<x2hA)B 

and 

®D2  “  “2^8 

Since  represents  a  zero-mean  white  noise  sequence, 
ECsiDiJ^jl  *  E ( S4oi  1 E C SlJ j  3  *  0  for  a11  i  *  j.  Further,  since 
the  present  state,  can  not  be  a  function  of  present  or 

future  inputs,  ^  is  independent  of  and  Kq2 •  By  the 
same  procedure  used  to  develop  the  discrete  covariance 
equation  given  by  equation  (19),  equation  (38)  will  be  used 
to  derive  the  discrete  covariance  equation  for  the  2nd 
order  Runge-Kutta  algorithm.  That  procedure  results  in 
PD(ti+1)  ■  5D(h,0)PD(ti)55(h,0) 

*  BdiQd1B5i  +  Bd2Qd2bS2  <39) 

where 

QDi  -  E(j£oi)fJlJ 

and 

Qq2  3  e<HD2e]52^ 

Mote  that  there  are  two  BqQqB^  terms  in  equation  (39). 

This  is  due  to  the  sampling  of  the  noise  generator  twice 
per  update,  once  per  derivative-function  evaluation.  If  it 
is  assumed  that  Qq  is  constant  over  the  update  interval, 
then  Qqi  ■  Qd2  *  Qd‘ 


Using  the  relationship  given  by  equation  (5),  it  is 
possible  to  define 

BqQqBjJ  =  BD1QDB5i  +  Bq2QqB($2 

From  the  definition  of  Bqj  and  Bq2  given  in  equation  (38), 
BDlQDBSl  =  ®^h^BQj)B^  +  a^c^a^h^tBQjjB^A1^  ♦  ABQqB  ) 

+  agh4(ABQDBTAT) 
and  Bq2QqB(£2  *  oc§h2BQQB^ 

Therefore,  for  the  2nd  order  Runge-Kutta  integrator, 

BqQqBJ  3  ( ccj  +  )h^BQjjB^  +  oc^o^a^h^CBQijB^A^  +  ABQqB^) 
+  o^h4(ABQDBTAT)  (40) 

Notice  that  BqQqB$  is  a  function  of  the  system  dynamics. 


represented  by  the  A  matrix.  Since  the  relationship  be¬ 
tween  Qq  and  Q  will  be  derived  using  BqQqB^  in  equation 
(22),  it  is  desirable  to  make  BqQqB^  independent  of  the 
dynamics  so  the  relationship  between  Qq  and  Q  will  be 
independent  of  the  dynamics.  To  accomplish  this,  the  same 
assumption  that  was  made  in  arriving  at  the  expression  on 
the  right  side  of  equation  (22)  will  be  made  again  to 
obtain  an  approximation  to  equation  (40)  that  is  indepen¬ 
dent  of  the  dynamics.  Specifically,  assume  that  h  is 
selected  reasonably  small,  such  that 


BDQDBg  s  <«i  +  «J)h2BQDBT 


For  oj  *  a2  =  0.5, 


BqQqBJ  “  h2BQDBT^i(«i)2  =  0.5h2BQDBT  (41) 

From  equations  (41)  and  (22),  it  is  seen  that  for  the  2nd 
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order  Runge-Kutta  integrator 

Qd  =  2Q/h 


(42) 


A  comparison  of  equation  (42)  to  equation  (33)  indicates 
that  the  relationship  between  Qq  and  Q  is  dependent  on  the 
integration  method. 

,Intsgg.gttff,r 

A  4th  order  Runge-Kutta  integrator  is  in  the  form  of 
*i+l  38  *i  +  <alk.l  ♦  «2^2  +  “3^3  +  (43) 

where 

tl  »  h£(2Li'  ih) ; 

)L2  m  hl<i.i  +  altl'(i  ♦  b1)h); 

It3  a  h£.(li  +  a2&-2'(i  +  b2)h); 

and  1^4  *  h£.(ii  +  a 3^3,  ( i  +  b3)h) 

Evaluating  the  j^'s  in  equation  (43)  for  the  system  given 


by  equation  (29)  yields 

Jil  «  h(Az,i  +  ®Hoi3  (44) 

lt.2  “  hI(I  +  a^hA)Aj.£  +  a^ABw^  +  Bj»q2  I  (45) 

ISL3  58  hid  +  a2hA  +  a2a jh3A2 ) Ajj.i  +  a2alh2^BHDl 

+  a2hABj£j)2  +  Bj£q3  3  (46) 

k>|  =  h[(I  +  a3hA  +  a3a2h2A2  +  a3a2aih3A3 ) A*! 

+  a3a2a1h3A3BH01  +  a3a2h2A2BsiD2 
+  a3hABst03  +  B&Q4 1  (47) 


where  the  Kd^'s  are  inputs  from  the  random  number  generator 
at  the  respective  evaluations  of  £(*,•).  Substituting 
equations  (44)  -  (47)  into  equation  (43)  results  in 


X.i+1  *  C(I  +  («1  +  0*2  +  c<3  +  OC4 )  hA 
+  +  a20c3  +  a3°^)h2A2 

+  (a2a^a3  +  a3a2°^4^^^^  +  ^a3a2al°^^4^4-IX-i 
+  Co^I  +  a^c^hA  +  a2aic«3h2A2 

+  a3a2aiO«4h2A23Bj£oi 

+  C«2l  +  a2oc3hA  +  a3a2c^h2A23Bj£02 

+  Ccigl  +  a  3  0(4  h  A  ]Bt£j)3  +  [o^IIBsq^  (48) 

A  popular  4th  order  Runge-Kutta  integrator  uses  = 

«4  *  1/6;  «2  =  cx3  =  1/3;  =  a3  =  =  b3  =  1/2;  and  a3  * 

b3  =  1.  Using  these  constants  to  evaluate  the  state  tran¬ 
sition  part  of  equation  (48)  results  in 

I  +  hi  +  0. 5(hA) 2  +  (hA) 3/6  +  (hA)4/24  =  5D(h,0)  (49) 

Note  that  equation  (49)  is  a  fourth-order  truncation  of 
3(h,0)  given  by  equation  (12).  Using  equation  (49),  equa¬ 
tion  (48)  can  be  rewritten  as 

JLi+1  *  $0 ( h , 0 ) li  +  BqiWqi  +  bD2*02  +  bD3S03  +  ®D4S04 
where 

®D1  “  C«il  ♦  ajo^hA  +  a2ai«3h2A2  +  a3a2al<x4h3A3^B' 
bD2  *  ^a2I  +  a2°*3hA  +  a3a2c<4h2A23B; 

Bq3  »  C0C3I  ♦  a3o^hA]B; 
and  Bq4  *  Co^IlB 

Equation  (50)  can  be  used  to  determine  the  second 
moment  function  for  this  algorithm.  A  number  of  statisti¬ 
cal  independence  relationships  will  simplify  the  math. 

Since  the  input  noise  is  a  zero-mean  white  noise  sequence, 
EtXQiK^j]  «  0  for  all  i  •  j .  Further,  the  present  state  is 


independent  of  present  and  future  inputs,  therefore, 

E[  iiiiSk  ]  =  0  for  k  *  1  to  4.  Using  these  statistical  rela¬ 
tionships,  the  discrete  covariance  matrix  for  the  4th  order 
Runge-Kutta  algorithm  can  be  derived  from  equation  (50)  by 
the  same  process  used  to  derive  equation  (19)  from  equation 
(17).  The  result  is 

pD(ti+l)  “  50(h,0)PD(ti)5j(h,0)  +  Bd1Qd1b{$i 

+  ®D2^D2®82  +  ®D3^D3®53  +  BD4^D4B54 
where  Qq^  =  E  ( SoitSiSk  1  for  k  =  1  to  4. 

Note  that  Pq^i+i)  has  four  Bq^Qq^B^  terms,  one  from  each 
sample  of  the  noise  generator  at  the  four  evaluations  of 
f(-,*). 

To  rewrite  equation  (51)  in  the  form  of  equation  (19), 
the  BojOd^bJi  terms  will  be  manipulated  into  one  equivalent 
BDQDBg  term.  To  do  this,  first  assume  that  Qq  is  constant 
over  the  update  interval,  leading  to  Oqj  *  Qq2  =  Qd3  =  Qd4 
=  Qq.  Now,  using  this  assumption  and  the  general  relation¬ 
ship  given  by  equation  (5),  define 

BD°DB!  “  BD1QdBSi  +  BD2^DB52  +  BD3^DB53  +  BD4Qdb84 
From  the  Bq^  definitions  given  in  equation  (50), 

BDiQDBSl  =  «fh2(BQDBT]  +  a^o^h3!  ABQDBT  +  BQDBTAT) 

+  a}a5h4(ABQDBTAT) 

+  a2a1a1a3h4IA2BQDBT  +  BQDBT(AT)2] 

+  a3a2ai«ia4hBCA3BQQBT  *  BQqBt(A^)3J 
♦  a  Ja2«2<x3h^  (  A2BQdBtAt  +  ABQDBT(AT)2] 

+  (a2a1oc3)2h6CA2BQDBT(AT)2l 
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+  a£a2a3cc2c<4h6[  A3BQqBtAt  +  ABQqBt(  At)  3  1 
+  aJa5a3«3c^h7(A3BQDBT(AT)2  +  A2BQDBT( AT) 3] 

+  (a3a2a1a4)2h8[A3BQDBT(AT)3] 

®D2®D®S2  "  °^h2CBQ[)B^]  t  a2®*2<*3h3 E ABQqB^  +  BQpB^A^] 

+  aJo«§h4tABQDBTAT} 

+  a3a2<x2cc4h4(A2BQDBT  +  BQDBT(AT)2] 

♦  a5a3o3o4h5EA2BQDBTAT  +  ABQDBT(AT)2J 
+  (a3a2o^)2h6EA2BQDBT(AT)2] 

BD3^DB53  =  °^h2EBQ|)B^]  +  a3cc3a^h3[  ABQqB^  +  BQqB^A-^] 

+  a§ajh4(ABQDBTAT] 

bD4Qdb54  =  «fr2[BQDBT] 

and  BqQqB^  *  BDlQDB5l  +  BD2®DB52  +  BD3QDBi$3  +  bD4Qdb54  (52) 
Notice  that  BqQqB§  is  a  function  of  the  system  dynamics  as 

it  was  with  the  2nd  order  Runge-Kutta  algorithm.  As  with 
the  2nd  order  algorithm,  it  is  desireable  to  make  BqQqB^ 
independent  of  the  dynamics  so  the  relationship  between  Qq 
and  Q  will  be  independent  of  the  dynamics.  To  accomplish 
this,  the  same  assumption  that  was  made  in  arriving  at  the 
expression  on  the  right  side  of  equation  (22)  will  be  made 
again  to  obtain  an  approximation  to  equation  (52)  that  is 
independent  of  the  dynamics.  Specifically,  assume  that  h 
is  selected  reasonably  small,  such  that 

BDQDBg  s  (of  ♦  ctj  +  o«5  ♦  aj)h2EBQDBT] 

For  ©n  “  a4  *  1/6  and  «  0*5  *  1/3 

bdQdb!  *  h2BQDBT^J«l)2  ■  h2BQDBT/3.6 
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With  this  approximation  to  BqQqB^  and  equation  (22),  the 
relationship  between  Qq  and  Q  for  this  4th  order  Runge- 
Xutta  integrator  is 

Qd  =  3. 60/h  (53) 


nth  Order  Runge-Kutta 

The  trend  established  in  the  above  analyses  can  be 
used  to  extend  the  results  to  an  nth  order  Runge-Kutta  in¬ 
tegrator.  Consider  that  the  coefficients  used  in  an  nth 
order  Runge-Kutta  integrator  are  selected  to  ensure  an 
nth  order  truncation  of  a  Taylor  series  expansion  of  the 
state  transition  matrix. [6]  Therefore,  the  integrator  will 
always  provide  a  5p ( h , 0 )  that  approximates  the  state  tran¬ 
sition  matrix,  5(h,0),  over  some  interval  h.  Based  on  this 
fact,  the  state  update  equation  for  the  system  will  be  in 
the  form  of 


*i+l  3  5D(h,0)i.i  +  j^BDfcMDfc 


(54) 


where 


BDk 


where 


p=n 


*m 


q  <  m 
q  =  m 


am'am+l *  *  *aq  >  m 
and  S{)(h,0)  is  an  nth  order  truncation  of  equation  (12). 

The  expression  in  equation  (54)  for  determining  Bq^  was 
derived  by  noticing  the  sequence  of  coefficients  in  the  BD 
terms  in  equation  (50),  equation  (38)  and  equation  (30). 


The  discrete  covariance  equation  for  a  system  solved 


by  an  nth  order  Runge-Kutta  will  be 

PD<ti+l>  =  5D<h'0)PD(ti>5?(h'0)  +  (55) 

where  QDJt  =  E C SLok.s5k. ^  *  Assuming  QDtc  =  QD  for  all  k  from  1 
to  n  and  assuming  that  h  is  small,  as  assumed  in  equation 
(21),  then,  for  an  nth  order  Runge-Kutta  integrator 

BdQdbS  *  h2BQDBT^i(ai)2  =  h2BQDBT/r  (56) 

where  r  »  l/c£  c£) 

From  equation  (56)  and  equation  (22),  the  relationship 
between  QD  and  Q  will  be 

Qd  *  TQ/h  (57) 

The  correction  factor,  r,  for  commonly  used  Runge-Kutta 
methods  is  given  in  Table  1. 


Table  1. 


Runge-Kutta  Correction  Factors 


Order 

coefficients 

Corr .  Factor 

2 

«1  -  -5  ,  a2  -  •  5 

3 

al  =  a3  “  <*2  =  2/3 

4 

=  1/6  ,  a 2  —  1/3 

ccg  =  1/3  ,  ocq  =  1/6 

4  (Gil) 

«1  -  1/6  ,  c«2  =  (2  -^2)/6 
CC3  =  (2  +v'2)/6  ,  04  =  1/6 

r  =  18/7 

=  2.57 

(5,6) 

ocj  =7/90  ,  «2  =0 
ag  =  32/90  ,014=  12/90 

0^  =  32/90  ,  ag  =  7/90 

r  =  3.54 

The  analysis  approach  taken  with  Runge-Kutta  methods 
becomes  difficult  to  manage  when  applied  to  other  integra¬ 
tion  methods.  Jury  (7]  offers  an  alternative  approach  that 
is  based  on  z-transforms  and  yields  a  steady-state  value  of 
the  autocorrelation  matrix.  Although  this  new  approach  is 
more  manageable  in  some  cases,  it  is  practically,  though 
not  mathematically,  restricted  to  scalar  equations.  To 
develop  and  understand  this  alternative  approach,  E(x2]  for 
a  2nd  order  Runge-Kutta  integrator  will  be  determined  by 
the  Jury  method. 

Jury  shows  that  if  a  discrete  scalar  system  driven  by 
discrete  white  noise  can  be  described  by 
M(z) 

X(z)  =  -  Wq( z )  *  G(z)Nq(z)  (58) 

L(z) 

where 

z”*  is  the  discrete  delay  operator; 

X(z)  is  the  z-transform  of  the  state  x^; 

H0(z)  is  the  z-transform  of  the  input  random 
sequence  with  variance  Qq; 

and  M(z)  and  Hz)  are  polynomials  in  z. 
then 

2  Qd  f*  ,  , 

Etx2]  =  -  0  G(z)G(z  1)z  (59) 

2ir  j  J  -it 

where  j  ■  /-I  * 

In  reference  [7],  Table  III  of  the  Appendix,  Jury 


provides  an  algorithm  for  evaluating  the  contour  integral 
in  equation  (59).  A  summary  of  the  algorithm  follows. 


Given  equation  (58)  such  that 


G(z ) 


M(z) 

L(z) 


'n  ♦  +  . . .  + 


n»jjz“  *  m^z* 


lgzn  +  ljZ 


n-1 


(60) 


♦  ...  +  l 


n 


then  equation  (59)  can  he  evaluated  to  yield  the  steady- 
state  value 

l«l| 


Elxf)  »  Qd- 


(61) 


lg|ft| 

where  |*|  denotes  the  matrix  determinant  operation.  The 
matrices  ft  and  %  are  square,  have  dimension  n+1  (n  is  the 
order  of  Hz)),  and  have  elements  formed  from  the  coeffi¬ 
cients  of  L(z)  and  M(z).  Specifically, 


ft  4 


lg  li 

11  l0+12 

12  13 


X2  *3 
11+13  12+l4 

10+14  1l4l5 


An 

*n-l 

ln-2 


(62) 


and  ft^  is  formed  by  replacing  the  first  column  of  ft  with 
the  vector 


To  apply  Jury's  algorithm  to  analyzing  the  2nd  order 
Runge-Kutta  integrator,  consider  a  scalar  linear  differen¬ 


tial  equation  of  the  form 

*(t)  =  Ax( t)  +  Bw(  t )  -  f(x,t)  (63) 

Applying  equation  (63)  to  the  Runge-Kutta  routine  given  by 
equations  (35)  -  (37)  yields  a  difference  equation  in  the 
form  of 

xk  “  *k-2  +  ®*Syk-i  +  0.5yk_2  (64) 

where 

xk-2  *  x(ih); 
xk  -  x( ( i+l)h) ; 
yk_2  “  hlAx(ih)  +  BwD1l; 
and  yk_x  -  h(A(yk_2  +  xk_2)  +  BwD2l 
where  wD1  and  wD2  are  defined  as  in  equation  (35).  Note 
that  k  increments  two  for  one  increment  in  ih.  This  is 
done  in  order  to  handle  the  two  sequential  inputs,  wD1  and 
wq2.  Taking  the  z-transform  of  equation  (64)  yields 
(0. 5z  +  0.5  *  0. 5hA)hBH(z) 

X(z)  »  - - -  (65) 

z2  -  5d 

where  3q  is  defined  above  equation  (38).  Note  that  equa¬ 
tion  (65)  is  in  the  form  of  equation  (58);  hence,  the  Jury 
algorithm  can  be  applied  directly.  For  this  case  mg  =  0; 
m^  *  0.5;  m2  ■  0.5(1  +  hA);  lg  ■  1;  lj  ■  0;  and  12  =  -5q. 
Evaluating  the  matrices  (i  and  &  using  these  coefficients 
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provides 


1-ID 


(0.5)2+(0.5)2(l+hA)2 
2(0. 5) 2( 1+hA)2 


1-5d 


Evaluating  the  determinates  of  ft  and  ft^  results  in 
|ft|  -  (1  -  5D)(1  -  Sg) 

and  |%|  -  (0.5)2(1  +  (1  +  hA)2l(l  -  5D) 

Therefore,  from  equation  (61) 


E(x{] 


(0.5  +  0 . 5 ( h A )  +  ( 0 . 5h A ) 2 ) h2B2QD 

1  -  5g 


(66) 


Manipulating  equation  (66)  yields 

E(xfl  -  jgEtxfl  +  (0.5  ♦  0. 5hA  +  (0.5hA)2)h2B2QD 
Coaparing  this  result  to  the  steady-state  value  of  the 
scalar  fora  of  equation  (19)  provides 

BgQo  *  l®-5  *  0 . 5hA  +  (0.5hA)2)h2B2QD 
which  agrees  with  the  scalar  form  equation  (40)  evaluated 
with  “  <*2  ■  0.5  and  a^  =  1. 


Adaas-Bashforth  methods  belong  to  a  class  of  integra¬ 
tion  algorithms  called  multistep  methods.  Multistep 
methods  have  the  form  [8] 


r  /v-  ov’o  v'?.  Z'~ 


to 


*i+l  •  ♦  ^ia3+l^(*-i-j'(i“j)h> 


(67) 


where  the  ctj's  ere  coefficients  dependent  on  the  particular 
algor i the  and  £(•,•)  is  given  in  equation  (29)  evaluated  at 
t  =  ( i-j )h.  Because  of  the  nature  of  eultistep  methods, 
the  analysis  approach  proposed  by  Jury  will  be  used  in 
determining  B§Qq. 

The  fourth-order  Adams -Bash forth  algorithm  for  the 
system  given  by  equation  (10)  is  (6) 

*i+l  “  Id  +  h(«1(Ax.i  +  ®SLi  >  +  «2<*Id-l  +  ZlLi-i* 

*  «3(Ali_2  +  B*d-2*  *  a4<Aj.i_3  ♦  BUj^))  (68) 
where  *  55/24  ;  aj  «  -59/24  ;  0*3  •  37/24  ;  *  -9/24  . 

Limiting  equation  (68)  to  the  scalar  case,  transforming  it 
to  the  z-domain  and  writing  it  in  the  form  of  equation  (58) 
yields 

hBCoiiZ3  +  aoz^  +  o^z  o41H(z) 

X(z)  -  - - -  (69) 

z*  -  (1  +  cxjhA )  z2  -  «2hAz2  -  o^hAz  -  04 h A 

where  X,  0,  A,  and  B  are  scalars.  The  application  of  equa¬ 
tion  (69)  to  the  algorithm  given  by  equations  (60)  -  (62) 
is  straight  forward,  but  quite  tedious,  because  it  requires 
solving  two  polynomial  matrices.  A  procedure  was  developed 
and  a  computer  program  was  written  to  help  facilitate  the 
the  computations.  That  procedure  and  program  are  provided 
in  Appendix  A.  The  results  of  the  Jury  Analysis  for  this 
problem  is 


„  h2B2QDCl  -  6. 34hA  +  18.9(hA)2  -  60.8(hA)3  +  ...3 

El*?]  -  - - 

1  - 

leading  to 

BgQo  -  h2B2QD Cl  -  6.34hA  +  18.9(hA)2  -  60.8(hA)3  +  ...3 

B$Qd  is  actually  a  polynomial  in  hA  of  infinite  order. 
The  first  fee  significant  terms  are  provided  above.  It 
will  be  assumed  that  |hA|  <  1,  thereby  making  it  possible 
to  approximate  the  B^Qq  by  a  finite  polynomial.  For  con¬ 
venience,  let 

bJQd  *  h2B2QDCl  -  6 . 34hA  +  18.9(hA)2  -  60.8(hA)33  (70) 

Using  equation  (70)  to  evaluate  the  left  side  of  the  scalar 
form  of  eqation  (22)  leads  to  the  relationship 

0D  -  Q/DKl  -  6. 34(hA)  ♦  18.9(hA)2  -  60.8(hA)3)3  (71) 


Notice  that  equation  (71)  has  a  strong  dependence  on  the 
system  dynamics.  However,  equation  (71)  will  reduce  to 
equation  (24)  as  hA  approaches  zero. 


IV.  ERROR  ANALYSES 


One  is  always  concerned  with  errors  in  performing 
analyses  o£  continuous  systems  via  simulation.  Digital 
simulation  introduces  a  number  of  unique  error  sources. 

The  effects  of  those  errors  are  usually  categorized  into 
two  areas;  one  pertaining  to  round-off  errors  (finite  word- 
length  errors)  and  one  pertaining  to  truncation  errors  due 
to  representing  the  system  models  by  truncated  infinite 
series.  With  respect  to  digital  simulations  of  physical 
systems,  round-off  errors  normally  begin  to  influence  the 
accuracy  of  the  results  when  the  integration  step  size,  h, 
is  decreased  to  the  point  that  the  changes  in  the  system 
dynamics  become  smaller  than  the  numerical  accuracy  of  the 
digital  computations.  On  the  other  hand,  truncation  errors 
typically  increase  as  the  step  size  becomes  larger.  Unac¬ 
ceptable  truncation  errors  occur  when  Sp(h,0)  no  longer 
accurately  approximates  5(h,0).  Usually,  there  is  a  region 
of  possible  values  of  h  in  which  neither  type  error  is 
significant.  For  efficiency  reasons,  the  system  analyst 
will  normally  select  the  largest  value  of  h  that  will  not 
introduce  significant  truncation  errors  and  use  a  machine 
that  provides  adequate  numerical  accuracy. 
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In  deterministic  analyses,  the  order  o£  the  integra¬ 
tion  method  is  usually  selected  to  minimize  truncation 
errors.  The  type  o£  nth  order  integration  method  is 
selected  to  maximize  efficiency  of  solution.  For  instance, 
it  is  generally  accepted  that  a  4th  order  Adams -Bash forth 
algorithm  is  more  efficient  than  a  4th  order  Runge-Kutta 
algorithm. (6]  In  stochastic  analyses,  these  same  problems 
exist  and  hence,  the  same  decisions  must  be  made,  but  an 
additional  error  source  is  introduced  due  to  the  need  to 
insure  statistical  accuracy.  In  Chapter  II,  it  was  shown 
that  the  first-order  statistics  of  the  noise  corrupted 
states  will  be  processed  in  the  same  way  as  deterministic 
states  but  higher  order  statistics  will  be  processed  in  a 
unique  manner.  The  second-order  statistics  will  be  proces¬ 
sed  by  a  solution  to  equation  (14)  and  higher  order  statis¬ 
tics  will  likewise  be  processed  by  unique  integral  equa¬ 
tions.  For  linear  systems  forced  by  Gaussian  distributed 
noise,  the  first  and  second-order  statistics  fully  describe 
the  statistical  behavior  of  the  system. (1)  However,  if  the 
noise  is  non-Gaussian  or  if  the  system  is  nonlinear,  higher 
order  statistics  will  have  to  be  calculated. 

Thus  far,  all  of  the  analyses  presented  in  this  paper 
with  regards  to  the  processing  of  the  second-order  statis¬ 
tics  have  been  based  on  the  assumption  that  h  would  be 
selected  small  enough  to  allow  one  to  approximate  the  state 
transition  matrix  by  the  identity  matrix.  In  practical 


terns,  this  means  that  the  step  size  required  to  accurately 
compute  the  second-order  statistics  must  be  much  smaller 
than  is  required  to  compute  the  first-order  statistics  as 
well  as  the  deterministic  states.  This  step  size  require¬ 
ment  will  translate  into  much  more  costly  analyses.  Fur¬ 
ther,  it  could  lead  to  the  introduction  of  round-off  er¬ 
rors.  Since  the  step  size  assumption  that  allowed  the  use 
of  equation  (21)  was  made  for  convenience,  it  will  be  use¬ 
ful  to  see  if  the  step  size  requirement  imposed  by  that 
assumption  can  be  relaxed.  To  perform  this  investigation, 
a  more  useful  form  of  equation  (20)  will  be  needed. 


Analysis  Approach 

If  one  restricts  the  problem  to  linear  systems  or  to 
well  behaved  (analytic)  nonlinear  systems  that  can  be  ac¬ 
curately  modeled  through  some  linearization  process,  then 
£(*,")  can  be  represented  by  the  infinite  series  given  in 
equation  (12).  Likewise,  the  integral  expression  in  equa¬ 
tion  (20)  can  be  represented  by  an  infinite  series.  To 
develop  the  series  representation  of  equation  (20),  let 


BQBt  =  A 


The  series  form  of  equation  (12)  is 

5<h''t,  *  E. — n— 


Using  equations  (73)  and  (74)  yields 

m  in  w  »  AiA(AT)^(h-T)i+j 

5(h,T)BQBT5T(h,T)  -  S  S :  - — - 

i»0  j»0  i  f  j  ! 


(73) 


(74) 


(75) 


Equation  (75)  is  a  polynomial  of  infinite  order  with  the 


independent  variable  T.  Since  A  and  A  are  independent  of 


T,  Equation  (75)  can  easily  be  integrated  with  respect  to  t 
and  evaluated  over  the  limits  0  to  h  as  in  equation  (20)  to 
yield 

f  5(h,t)A5T(h,T)dT  -  S  S  [(h-T)i+*dT 

.  -S?  JS  AiA(AT)3(h)i+*  _ 

*  h  S  SI  -  (76) 

i-0  j =0  ( i+j+1) i  !  j ! 


Equation  (76)  is  a  quadratic  polynomial  in  hA  of  in¬ 
finite  order.  If  h  is  selected  reasonably  small  such  that 
h*|Xmax|  <  1*  where  |Xmar |  is  the  magnitude  of  the  largest 
eigenvalue  of  A,  then  equation  (76)  can  be  reasonably  ap¬ 
proximated  by  a  finite  series  of  order  n.  A  truncated 
expression  of  equation  (76)  can  be  used  to  determine  the 
order  of  magnitude  of  the  local  error  in  the  discrete  co- 
variance  computations.  The  appropriate  order  of  truncation 
will  be  a  function  of  the  integration  method  and  its  selec¬ 
tion  will  be  based  on  the  order  of  the  BqQqB§  expression 
determined  for  the  integration  method.  Thus,  at  this 
point,  the  analysis  becomes  integration  method  dependent. 

Integrator  Analyses 

In  the  analyses  to  follow,  equation  (76)  will  be 
truncated  at  various  orders  in  h.  For  notational 
brevity,  an  nth  order  truncation  of  equation  (76) 
means  that  i  will  step  from  0  to  n  and  j  will  step 
from  0  to  n,  however  all  terms  resulting  from  the 
condition  ((i+j)  >  n]  will  be  neglected. 


(77) 


For  the  Euler  integrator,  from  equation  (32), 

BDQDb5  =  h2BQDBT 
Equating  the  right  side  of  equation  (77)  to  a  first-order 
truncation  of  equation  (76),  yields 

h2BQDBT  -  h(BQBT  +  0.5h(ABQBT  +  BQBTAT))  (78) 

If  Qq  =  Q/h,  as  proposed  by  Griffith[4),  then  equation 
(78)  will  be  valid  if  the  first-order  terms  in  h  on  the 
right  side  of  equation  (78)  are  much  smaller  than  the 
zeroth-order  term.  However,  if  it  is  assumed  (as  Griffith 
did)  that  the  state  transition  matrix  can  be  approximated 
by  the  identity  matrix  (a  zeroth-order  truncation  of  equa¬ 
tion  (12))  then  equation  (24)  will  be  a  valid  relationship 
that  can  be  used  to  satisfy  equation  (78).  Making  that 
assumption  translates  into  a  smaller  step  size  requirement. 
Therefore,  it  is  concluded  that  the  local  error  that 
results  in  using  equation  (24)  to  relate  Qq  to  Q  for  the 
Euler  integrator  is  of  order  h  as  compared  to  a  local 
integration  error  of  order  h2.[6] 

2nd  Order  Runoe-Kutta 

Hith  =  o*2  =  0.5,  equation  (40)  evaluates  to 

BDQDB{  -  0.5h2(BQDBT  +  0.5h(ABQDBT  +  BQDBTAT ) 

+  h2(ABQDBTAT/2) )  (79) 

The  right  side  of  equation  (76)  truncated  at  2nd  order  is 
h(BQBT  +  0.5h(ABQBT  +  BQBTAT)  +  h2(ABQBTAT/3 
+  (A2BQBt  ♦  BQBT(AT)2)/6) ) 


(80) 


By  direct  inspection,  it  is  seen  that  the  right  side  o£ 
equation  (79)  will  approximate  expression  (80)  with  local 
error  on  the  order  of  h2,  if  Op  =  2Q/h.  This  compares  to  a 
local  integration  error  on  the  order  of  h3.(6] 


4th  Order  Runoe-Kutta 

For  notational  convenience  let 

A  -  BQBt 
and  /^)  »  BQdBt 

For  the  most  commonly  used  4th  order  Runge-Kutta  integrator 
with  ai  ■  <*4  *  1/6;  <*2  ■  o«3  «  1/3;  a^  =  a2  =  0.5;  and  a2  = 
1,  equation  (52)  yields 


BdqDb5  *  +  +  A|)AT3h  +  ~--“-CA2Ao  +  /^(A7)2^2 

+  3^33CAA°Al:"’2  *  «7“3Ad  *  /'0uT,3:lh3 

t  -i-CA2/\j,AT  ♦  AAo<AT)2Jh3  *  -i-CA2Ao<AT)2]h4 

♦  “~“CA3/\qAt  ♦  A/\o(AT)33h4 

♦  r^~CA3/\0(AT)2  +  A2/^(AT)33h5 
80  • 

♦  ^-CA3/y)(AT)3Dh6  ]  (81) 

160  J 

The  right  side  of  equation  (76)  truncated  at  sixth-order  is 


h£  A  +  . 5  CAA  +  AAT3h  +  -£-CA2A  +  A(AT)23h2 

♦  -^-CAAAT3h2  +  — t~CA3A  +  A(AT)33h3 

3  24. 


1. 

3 
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♦  —  CA2AAT  +  AA(AT)23h3  +  r^-CA2A(  AT)  23h4 

8  20. 

♦  r^“CA3AAT  ♦  AA(AT)33h4  +  r~CA4A  +  A(AT)43h4 

30.  120 

♦  :TTca5A  +  A(AT)53h5  +  CA3A(At)2  ♦  A2A(AT)33h5 

720  72. 

+  ^-£A4AAT  ♦  AA(AT)43h5  +  ^jCA3AD(AT)3Dh6 

♦  tA6A  +  A(AT)63h6  +  ~—CA4A(At)2  +  A2A(AT)43h6 

5040.  336 


^■CA5AAT  *  AA(AT)53h6| 


(82) 


A  comparison  of  the  right  side  of  equation  (81)  to  expres¬ 
sion  (82)  reveals  that  equation  (81)  will  approximate  ex¬ 
pression  (82),  with  local  error  on  the  order  of  h2,  if 

3. 6Q 

°» — r~ 

Note  that  a  4th  order  Runge-Kutta  has  an  error  of  the 
same  order  in  h  as  a  2nd  order  Runge-Kutta.  Based  on  this 
observation,  one  may  ask  if  this  means  that  a  2nd  order 
Runge-Kutta  should  be  preferred  over  a  4th  order  Runge- 
Kutta  for  efficiency  versus  accuracy  reasons.  Certainly 
not,  for  two  reasons.  First,  the  4th  order  Runge-Kutta  is 
still  preferred  for  the  deterministic  state  equations  and 
for  the  first-order  statistics  of  the  non-deterministic 
states.  Second,  although  the  two  integrators  have  covari¬ 
ance  errors  on  the  same  order  in  h,  the  absolute  error  of 
the  4th  order  Runge-Kutta  will  be  smaller  than  that  of  the 
2nd  order  Runge-Kutta.  To  see  this,  consider  that  the  abso 
lute  error  in  the  4th  order  integrator  is  a  term  by  term 
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difference  between  equation  (81)  and  equation  (76).  Like¬ 
wise,  the  absolute  error  in  the  2nd  order  integrator  is  a 
tern  by  term  difference  between  equation  (79)  and  equation 
(76).  For  comparison  purposes,  expression  (81)  contains 
all  the  useful  information  of  equation  (76)  because  both 
methods  are  equally  in  error  above  order  six.  Since  the 
second-order  integrator  contains  no  information  above  order 
two  and  the  4th  order  integrator  only  deviates  by  a  frac¬ 
tional  amount  in  each  common  term  through  order  six,  it  is 
concluded  that  the  4th  order  method  has  a  smaller  absolute 
error. 

flth.  .Order  Adama -Bash forth  Integrator. 

From  equation  (70),  the  scalar  case  of  the  4th  order 
Adams-Bashforth  integrator  provides 

BgQD  *  h2B2QD Cl  -  6. 34hA  +  18.9(hA)2  -  60.9(hA)3l  (83) 

A  third-order  truncation  of  the  right  side  of  equation  (76) 
for  the  scalar  case  yields 

hB2QCl  +  hA  +  2(hA)2/3  +  (hA)3/33  (84) 

Since  equation  (83)  and  expression  (84)  are  scalar  func¬ 
tions,  the  relationship  between  Qq  and  Q  can  be  determined 
directly.  Doing  so,  yields 

QC1  ♦  hA  +  2(hA) 2/3  +  (hA)3/33 

Qd  »  - - - - -  (85) 

h Cl  -  6 . 34hA  +  18.9(hA)2  -  60.9(hA)3l 

Casual  inspection  of  equation  (85)  shows  that  equation  (24) 

is  only  valid  for  hA  «  1  when  applied  to  the  Adams- 

Bashforth  integrator. 
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It  will  be  beneficial  to  examine  the  above  results  for 
the  scalar  case.  Restricting  the  equations  to  scalar  sys¬ 
tems  will  allow  the  results  to  be  demonstrated  graphically 
by  defining  an  error  equation  and  plotting  it  as  a  function 
of  the  parameter,  hA.  The  parameter,  hA,  is  the  product  of 
the  step  size,  h,  with  the  dynamic  coefficient,  A.  For  the 
scalar  case,  A  is  the  reciprocal  of  the  system  time-con¬ 
stant.  The  step  size  is  selected  such  that  the  ratio  of  h 
to  the  system  time-constant  is  less  than  one,  typically 
one-tenth.  Therefore,  hA  is  typically  0.1  or  less. 

A  scalar  case  of  equation  (20)  shows  that,  ideally, 
for  each  integrator  considered,  the  normalization  function 
should  result  in  the  equality 

bBqD  ■  [  52(h,DB2QdT 

*0 

fh 

=  B2Q  €2A(h-T)dt 

<0 

-  B2Q(€2hA  -  1)/2A  (86) 

In  each  of  the  cases  presented  in  this  chapter, 

Bg  =  h2B2f (hA)/T 

where  f(hA)  was  a  polynomial  in  hA,  and  r  was  a  constant 
used  to  make  the  zeroth  coefficient  of  f(hA)  equal  to  1. 
Further,  in  order  to  make  the  relationship  between  Qq  and  Q 
independent  of  the  system  dynamics,  it  was  proposed  that 
for  each  of  the  cases,  with  the  possible  exception  of  the 


Adans -Bash forth  integrator. 


Therefore, 


Qd  =  TQ/h 


bJ0D  -  B2Qhf(hA)  (87) 

Substituting  equation  (87)  into  equation  (86),  provides  the 
ideal  relationship 

2hAf(hA)  -  €2hA  -  1 

In  reality,  this  relationship  will  not  be  exact  because 
f(hA)  *  ( f 2hA  _  i)/2hA 


Based  on  this  realization,  an  error  function  can  be 
defined  and  evaluated  for  each  integrator.  Let 

V(hA)  =  2hAf(hA)  -  €2hA  +  1  (88) 

where  V(hA)  is  the  error  function.  Note  that  V(hA)  *  0  for 
hA  *  I.  This  function  can  be  plotted  versus  hA  for  each  of 
the  f(hA)  functions  derived  in  these  analyses.  A  number  of 
practical  assumptions  will  be  made  in  plotting  V(hA)  for 
each  of  the  integrators.  First,  it  will  be  assumed  that 
the  dynamic  system  being  simulated  is  stable;  thus,  A  <  0. 
Next,  the  range  of  hA  will  be  restricted  to  values  that 
would  likely  be  considered  for  use  in  a  simulation;  specif¬ 
ically,  0  s:  | hA |  £  0.2. 

Euler  Integrator 

Prom  equation  (77),  for  the  Euler  integrator,  it  is  seen 
that  f(hA)  *  1.  Figure  2  provides  a  plot  of  equation  (88) 
evaluated  with  f(hA)  ■  1.  Note,  that  using  equation  (24) 
to  set  the  random  number  generator's  covariance  would 


For  the  2nd  order  Runge-Kutta  integrator,  equation 


(79)  provides 

f (hA)  *  1.0  +  hA  +  0. 5(hA) 2 

Likewise,  for  the  4th  order  Runge-Kutta  integrator,  equa¬ 
tion  (81)  provides 

f (hA)  -  1.0  +  hA  +  (3/5) (hA)2  +  (l/4)(hA)3 

♦  (l/10)(hA)4  +  (l/40)(hA)5  +  (l/160)(hA)6 


Figure  3  provides  a  comparison  plot  of  V(hA)  evaluated  with 
these  two  f(hA)  functions.  Note  that,  as  predicted  in  the 
error  analyses,  the  order  of  magnitude  of  the  error  is  the 


same  for  the  two  methods  and  the  absolute  error  of  the  4th 
order  Runge-Xutta  Integrator  is  smaller  than  that  of  the 
2nd  order  Runge-Kutta  integrator.  However,  both  integra¬ 
tors  provide  accurate  results  for  the  second-order  statis¬ 
tics  when  the  proper  statistical  relationship  is  used. 
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There  are  several  choices  for  the  Qq  to  Q  normalizing 
function.  Consider  two  choices.  First,  consider  the  one 
proposed  by  Griffith  where  Qq  =  Q/h.  The  advantage  of  this 
choice  is  that  the  normalizing  function  is  independent  of 
the  system  dynamics.  With  this  normalizing  function, 
f (hA)  *  1.0  -  6 . 34hA  ♦  18.9(hA)2  -  60.9(hA)3 
The  second  choice  would  be  to  use 

Qq  -  Q/(h(1.0  -  6. 34hA  +  18.9(hA)2  -  60.9(hA)3)) 

and 

f (hA)  *  1.0 

The  obvious  disadvantage  of  this  choice  is  that  the  Qq  to  Q 
normalizing  function  is  dependent  on  the  system  dynamics. 
However,  for  analysis  purposes,  this  choice  must  be  consi¬ 
dered.  A  comparison  plot  for  these  two  f(hA)  functions  is 
given  in  Figure  4.  Note  from  Figure  4,  that  the  use  of 
Qq  =  Q/h  for  the  Adams -Bash forth  integrator  could  result  in 
errors  on  the  order  of  100  percent.  Further,  note  that 
even  when  using  the  complex  normalizing  function,  the  er¬ 
rors  are  the  same  as  would  be  obtained  from  using  an  Euler 
integrator.  Besides  this  sensitivity  to  step  size,  there 
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are  other  problems  with  the  use  of  the  Adams-Bashforth 
integrator  for  stochastic  analyses.  These  other  problems 
will  be  demonstrated  in  Chapter  V. 


V.  NUMERICAL  AND  FURTHER  ANALYSES 


To  demonstrate  the  significance  of  the  findings  in 
Chapters  III  and  IV,  a  first-order  linear  differential  sto¬ 
chastic  equation  will  be  analyzed  numerically.  Limiting 
the  example  to  a  linear  scalar  system  will  allow  the  covar¬ 
iance  to  be  conveniently  expressed  analytically  and  the 
numerical  results  compared  directly.  For  the  example,  the 
continuous  system  is  described  by 

x(t)  «  -2x(t)  +  w ( t )  (89) 

where  w(t)  is  a  normally  distributed  white  random  process 
with  a  spectral  density  equal  to  e2.  For  the  example  let 
c2  «  .01  . 

Equation  (14)  will  be  used  to  analytically  determine 
P(t).  At  the  outset,  assume  that  tg  *  0  and  P(tg)  =  0 . 

With  these  assumptions,  the  scalar  form  of  equation  (14) 
applied  to  this  system  is 

P(t)  =  [  52(t,T)B2Q(T)dt 

*2 

where  S(t,T)  »  c"2**”*1*  ;  B  =  1;  and  Q(D  ■  r2.  Evaluating 
the  right  side  of  the  covariance  equation  provides 

P(t)  «  v2(l  -  €"4t)/4  (90) 

To  numerically  analyze  this  system,  a  Monte  Carlo 
simulation  was  written  in  FORTRAN  77  and  executed  using 


single  precision  arithmetic  on  a  VAX  780  computer.  One 
thousand  executions  were  made  per  Monte  Carlo  simulation. 
The  simulation  was  written  so  the  integrator  would  be  a 
module  interfaced  with  the  rest  of  the  simulation  through 
parameter  I/O.  Each  integrator  of  interest  was  then  coded 
as  a  module  and  incorporated  in  the  simulation  as  needed. 
Thus,  the  only  variable  from  one  simulation  program  to 
another  was  the  explicit  integrator  under  investigation. 

A  random  number  generator  proposed  by  Marse  and  Roberts  (9) 
was  used  to  generate  a  uniform  distributed  random  sequence. 
The  uniformly  distributed  sequence  was  then  shaped  into  an 
approximate  Gaussian  distributed  sequence  by  direct  appli¬ 
cation  of  the  Central  Limit  Theorem  involving  the  summation 
of  twelve  samples  from  the  random  number  generator. 11,21 
The  random  number  generator  routine  is  given  in  Appendix  B. 

Method  Results 

Figures  5  through  10  show  the  results  from  the  Monte 
Carlo  analyses  for  various  fixed-step  integrators.  The 
integrators  were  selected  to  exercise  and  test  the  analyti¬ 
cal  result  given  by  equation  (57)  for  a  number  of  frequent¬ 
ly  used  Runge-Kutta  algorithms.  Table  2  provides  a  summary 
of  the  methods  tested  as  well  as  the  step  size  and  normali¬ 
zation  function  used  in  each  test. 


Table  2 


Summary  of  Figures  Providing  Fixed-Step  Results 


INTEGRATOR 

METHOD 

STEP  SIZE 
( seconds ) 

NORMALIZING 

FUNCTION 

5 

EULER 

0.005 

QD  *  Q/h 

6 

2nd  Ord  R-K 

0.005 

Qd  =  2Q/h 

7 

4th  Ord  R-K 

0.005 

Qd  =  3 . 6Q/h 

8 

4th  Ord  R-K 

0.100 

Qq  =  3 . 6Q/h 

9 

4th  Order 
R-K -Gil 

0.050 

Qd  =  18Q/7h 

10 

R-K  (5,6) 

0.050 

QD  =  3 . 54Q/h 

The  step  size  used  to  generate  the  data  in  Figures  5, 
6,  and  7  was  selected  to  enforce  the  original  assumption 
that  the  state  transition  matrix  would  approximate  the 
identity  matrix.  At  a  step  size  of  0.005  seconds,  the 
state  transition  matrix  for  the  system  described  by  equa¬ 
tion  (89)  is  approximately  equal  to  0.99  over  one  step.  The 
plots  clearly  demonstrate  the  validity  of  equation  (57). 

The  results  also  show  that  the  use  of  Qq  -  Q/h  for  any 
Runge-Kutta  method  higher  than  first-order  will  introduce 
significant,  and  unnecessary  errors  into  the  simulation. 

The  error  analyses  in  Chapter  IV  predicted  that  the 
normalization  function  for  the  Runge-Kutta  methods  would  be 
fairly  insensitive  to  a  relaxation  of  the  step  size  re¬ 
strictions  imposed  by  the  assumption  that  the  state  trans¬ 
ition  matrix  approximate  the  identity  matrix.  To  test  that 
prediction,  the  4th  order  Runge-Kutta  integrator -based 
simulation  was  re-run  at  a  large  step  size  of  0.1  seconds. 
Figure  8  shows  the  results  from  this  test.  It  should  be 
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noted,  that  although  this  seems  to  be  a  large  step  size, 
for  this  example, 

|2d(0.1,0)  -  5(0. 1,0) J  <  3*10“® 
where  5(*,«)  is  given  by  equation  (12)  and  5q(*,*)  is  given 
by  equation  (49).  As  seen  by  the  plot,  the  normalization 
function  still  provides  accurate  numerical  results. 

Figure  9  shows  the  results  for  a  Runge-Kutta-Gil  meth¬ 
od.  [5]  The  Runge-Kutta-Gil  method  is  a  4th  order  method 
that  uses  a  non-convent ional  algorithm  for  calculating  the 
intermediate  update  vectors.  Gil  derived  his  procedure  to 
minimize  finite  word-length  errors.  The  algorithm  was 
tested  to  insure  that  equation  (57)  would  hold  for  this 
nonconvent ional  method.  Clearly  it  does. 

Figure  10  shows  the  results  for  another  non-conven- 
tional  method.  This  method  is  called  a  Runge-Kutta  (5,6). 
[5]  It  is  a  sixth-order  method  and  thus  makes  six  evalua¬ 
tions  of  the  derivative  functions  during  each  update.  How¬ 
ever,  it  only  explicitly  uses  five  of  the  intermediate  up¬ 
date  vectors  to  propagate  the  states  over  the  interval . 

This  algorithm,  as  well  as  the  Runge-Kutta  (7,8)  method,  is 
often  used  in  situations  needing  high  accuracy,  such  as  the 
solution  to  orbital  equations  of  motion.  The  correction 
factor,  T,  for  the  Runge-Kutta  (5,6)  integrator  is  given  in 
Chapter  III,  Table  1.  Figure  10  shows,  once  again,  that 
equation  (57)  is  a  valid  general  solution  to  the  problem. 


The  Adams -Bash forth  integrator  selected  for  evaluation 
is  a  fourth-order  method  and  is  given  by  equation  (68). 

The  integrator  was  first  tested  at  an  integration  step  of 
0.005  seconds.  Equation  (71)  was  used  to  determine  the 
correct  value  of  the  discrete  covariance.  The  results  of 
this  numerical  analysis  are  provided  in  Figure  11.  Figure 
11  shows  a  comparison  between  the  calculated  state  covari¬ 
ance  when  equation  (71)  was  used,  the  calculated  state  co- 
variance  when  equation  (24)  was  used,  as  suggested  by 
Griffith,  and  the  analytical  solution  to  the  state  covari¬ 
ance  given  by  equation  (90).  Figure  11  demonstrates  that 
equation  (71)  provides  greater  accuracy  over  equation  (24). 
Examination  of  equation  (71)  shows  that  as  the  step  size 
approaches  zero,  equation  (71)  will  approach  equation  (24). 
Figure  12  demonstrates  this  relationship  and  provides  val¬ 
idity  to  Griffith's  numerical  findings.  For  the  case  pre¬ 
sented  in  Figure  11,  the  step  size  was  selected  to  be  0.001 
seconds . 

An  interesting  problem  occurred  with  the  Adams-Bash- 
forth  integrator.  This  problem  occurs  only  in  stochastic 
analyses.  Recall  that  the  Adams -Bashf or th  method,  as  well 
as  all  multistep  methods,  are  not  self-starting.  To  get  an 
nth  order  multistep  integrator  started,  computation  of  the 
derivative  function  over  the  first  (n-1)  steps  must  be 
accomplished  independent  of  the  multistep  integrator. 


Figure  11.  4th  Order  fldaas*6ashforth  Analysis  (h  -  0.005  seconds) 
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Often  a  4th  order  Runge-Kutta  method  is  used  to  start  a  4th 
order  Adams -Bash forth  integrator.  The  order  of  the  two 
methods  is  chosen  to  be  the  same  so  that  the  numerical  ac¬ 
curacy  is  of  the  same  order  of  magnitude.  The  problem  that 
occurs  in  stochastic  analyses  using  the  Adams -Bash forth 
method  is  due  to  the  starting  procedure. 

In  order  for  the  starting  integrator  to  provide  accu¬ 
rate  statistics  over  the  starting  interval,  the  covariance 
of  the  random  number  generator  must  be  set  with  respect  to 
the  starting  integrator.  For  instance,  if  a  4th  order 
Runge-Kutta  algorithm  is  used  to  start  the  multistep  inte¬ 
grator,  Qq  must  be  related  to  Q  by  Qq  =  3.6Q/h  over  the 
starting  interval.  Once  the  four  evaluations  have  occur¬ 
red,  control  is  turned  over  the  Adams -Bash forth  integrator 
and  Qq  is  then  set  in  accordance  with  equation  (71).  The 
problem  occurs  during  the  first  four  updates  made  by  the 
Adams -Bash forth  integrator.  The  noise  terms  used  in  the 
evaluation  of  the  derivative  functions  during  the  starting 
interval  are  mismatched  to  the  noise  terms  needed  by  the 
Adams -Bash forth  for  accurate  propagation  of  the  states. 

This  mismatch  causes  a  large  transient  spike  and  introduces 
an  error  in  the  state  calculation  that  propagates  in  time. 
Figure  13  illustrates  this  phenomenon.  This  problem 
appears  only  in  the  calculation  of  the  stochastic  states. 
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To  remedy  the  problem,  a  number  o£  potential  solutions 
was  tried  resulting  in  varying  degrees  of  success.  The 
necessary  requirements  for  an  adequate  solution  were  found 
to  be  that  the  starting  procedure  must  provide  accurate 
updates  at  the  selected  step  size  and  the  noise  normaliza¬ 
tion  function  used  by  the  starting  integrator  must  be  ap¬ 
proximately  the  same  as  the  one  used  by  the  Adams -Bash forth 
The  method  that  finally  provided  good  results  was  a  method 
involving  the  use  of  a  number  of  different  integration 
methods.  At  the  first  update  an  Euler  method  was  used, 
followed  by  a  2nd  order  Adams-Bashforth  integrator  at  the 
second  update,  followed  by  a  3rd  order  Adams-Bashforth 
integrator  at  the  third  update  and  finally  the  4th  order 
Adams-Bashforth  integrator  at  the  4th  and  subsequent  up¬ 
dates.  The  data  used  to  construct  Figures  11  and  12  were 
determined  with  the  4th  order  Adams-Bashforth  integrator 
started  by  this  method. 

The  error  analyses  in  Chapter  IV  predicted  that  the 
Adams-Bashforth  integrator  would  be  sensitive  to  a  viola¬ 
tion  of  the  assumption  that  the  state  transition  matrix 
would  approximate  the  identity  matrix  over  the  integration 
interval.  To  test  this  prediction,  the  Adams-Bashforth 
simulation  was  run  at  a  step  size  of  0.025  seconds.  At 
this  step  size  3(h,0)  a  0.95  .  This  step  size  was  chosen 
because  the  total  execution  time  realized  for  the  simula¬ 
tion  was  138  seconds  which  is  slightly  longer  than  that 
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realized  running  the  4th  order  Runge-Kutta  simulation  (126 
seconds)  with  a  step  size  o£  0.1  seconds.  Figure  14  shows 
the  results  from  this  test  and  verifies  the  prediction.  A 
Comparison  of  Figure  14  to  Figure  8  and  considering  equiva¬ 
lent  execution  times,  indicates  that  there  is  no  benefit  in 
using  an  Adams -Bash forth  integrator  for  stochastic 
simulations. 


Figure  13.  4th  Order  Rdais-Bashforth  with  Runge-Kutta  Starter 
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Figure  14.  fldaas-Bashforth  Analysis  (h  -  0.025  seconds) 
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Thus  far,  all  of  the  analyses  have  dealt  with  conven¬ 
tional  integrators.  A  number  of  more  complex  integration 
methods  are  frequently  used  in  special  situations.  These 
more  complex  methods  are  usually  selected  in  order  to 
obtain  highly  accurate  solutions  at  moderate  computational 
costs.  Two  types  of  "advanced"  algorithms  that  are  of  in¬ 
terest  to  this  study  are  variable-step-size  methods  and 
predictor -corrector  methods. 


-1  -1  •»  ,  C-M  jW<  P 


Variable-step-size  methods  are  useful  for  simulating 
systems  that  contain  time-varying  dynamics.  In  any 


simulation,  the  step  size  should  be  selected  to  handle  the 
fastest  dynamics  in  the  system.  However,  if  the  dynamics 
are  fast  over  one  interval  and  then  slow  over  another 


interval,  a  step  size  selected  for  the  fast  interval  will 
cause  wasted  computations  during  the  slow  interval.  Varia¬ 
ble-step-size  methods  address  this  problem  by  computing 
some  function  that  senses  the  accuracy  of  solution  and 
adaptively  adjusts  the  step  size  so  the  errors  remain 
bounded  below  some  defined  tolerance  level.  Application  of 
variable-step-size  methods  to  stochastic  systems  presents 
some  unique  problems  not  encountered  with  deterministic 
systems . 

As  established  in  this  dissertation,  the  magnitude  of 
the  noise  added  to  the  system  states  is  a  function  of  the 
step  size.  This  fact  translates  into  effectively  discount¬ 
ing  all  variable-step-size  methods  based  on  multistep  in¬ 
tegration  methods  because  a  step  size  adjustment  would 
require  a  total  restart  of  the  integrator  using  new  noise 
terms.  Therefore,  at  the  outset,  the  search  for  candidate 
variable-step-size  methods  for  this  analysis  will  be  re¬ 
stricted  to  Runge-Kutta  methods.  As  seen  in  the  analyses 
presented  in  Chapter  III,  each  fixed-step  integrator  has  a 
unique  normalizing  function  for  setting  the  random  number 
generator's  covariance.  This  fact  makes  the  application  of 
methods  similar  to  the  Runge-Kutta-Fehlburg  integrator  dif¬ 
ficult  to  manage. 


The  Runge-Kutta-Fehlburg  integrator  is  an  efficient 
algorithm  that  calculates  the  state  update  by  solving  a 
fourth-order  Runge-Kutta  integrator  and  a  fifth-order 
Runge-Xutta  integrator  in  parallel.  The  fourth-order  solu¬ 
tion  is  used  to  propagate  the  states  and  the  fifth-order 
solution  is  used  to  predict  the  error  of  the  fourth-order 
solution.  The  efficiency  of  the  Fehlburg  method  is  found 
in  the  intermediate  derivative  evaluations.  A  fourth-order 
Runge-Kutta  method  requires  four  function  evaluations  and  a 
fifth-order  Runge-Kutta  method  requires  six  function  evalu¬ 
ations.  Therefore,  a  brute  force  implementation  o ?  a 
fourth  and  a  fifth-order  method  would  require  ten  function 
evaluations.  Fehlburg  derived  a  set  of  coefficients  that 
would  allow  the  fourth-order  update  formula  to  use  the  same 
derivative  function  evaluations  as  the  fifth-order  update; 
thus  the  entire  process  only  requires  six  function  evalua¬ 
tions  as  compared  to  ten  by  the  brute  force  method. (6] 

The  difficulty  with  applying  Fehlburg's  algorithm  to 
stochastic  systems  is  due  to  the  need  to  use  one  noise 
covariance  for  the  fourth-order  update  and  another  for  the 
fifth-order  update.  This  means  that  the  two  update  formu¬ 
las  can  no  longer  share  the  derivative  function  evalua¬ 
tions.  The  requirement  of  ten  function  evaluations  per 
step  makes  the  computational  costs  prohibitively  large. 

Any  variable-step-size  method  that  uses  two  update  methods 


of  different  order  to  calculate  the  error  function  would 


introduce  this  same  problem. 

A  variable-step-size  method  that  does  not  have  any  of 
the  above  problems  is  one  presented  in  reference  5.  This 


method  was  developed  by  Colatz  and  is  based  entirely  on  a 
4th  order  Runge-Kutta  integrator.  The  error  formula  that 
Colatz  derived  is  a  function  of  the  intermediate  function 
evaluations.  Specifically,  the  error  function,  €,  is 


k2  “  K3 
*1  -  *2 

wfc-*re  the  Xj/s  are  given  in  equation  (43). 


(91) 

Colatz  showed 


that  the  error  function  given  by  equation  (91)  should  be 


less  than  a  few  hundredths  in  order  to  insure  an  accurate 


solution.  The  only  difficulty  in  applying  Colatz 's  formula 
was  the  situation  where  all  the  X^'s  are  nearly  zero.  This 
situation  occurs  in  steady-state.  Explicit  implementation 
of  equation  (91)  in  computer  code  will  cause  numerical  er¬ 


rors  in  this  situation  due  to  a  division  by  small  numbers. 
To  circumvent  that  problem,  equation  (91)  was  implemented 
in  the  following  form. 

|X2  -  X3|  -  CJXj  -  X2|  *  0  (92) 

where  €  was  set  to  0.01  .  When  equation  (92)  was  used  to 
predict  the  maximum  step  size  for  a  deterministic  system 
with  the  same  dynamics  as  in  equation  (92),  it  predicted 
that  h  <  0.1  seconds.  This  is  certainly  a  reasonable  fin¬ 
ding.  However,  when  equation  (92)  was  applied  to  the  sto- 
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chastic  system  given  by  equation  (89),  the  results  were 
clearly  erroneous.  Figure  15  is  a  plot  o£  the  left  side  of 
equation  (92)  versus  time  for  the  stochastic  system  solved 
at  a  step  size  of  0.005  seconds.  As  seen  from  equation 
(92),  the  evaluated  expression  must  be  less  than  or  equal 
to  zero  in  order  to  meet  the  accuracy  condition.  Figure  15 
shoes  that  this  condition  is  never  met  indicating  a  need  to 
reduce  the  step  size.  Although  other  analyses  have  shown 
that  this  step  size  is  well  within  acceptable  limits,  the 
step  size  was  reduced  an  order  of  magnitude  to  an 
h  *  0.0005  seconds.  Figure  16  shows  the  computed  error 
function  at  this  step  size.  Clearly,  from  Figures  15  and 
16,  the  addition  of  noise  to  the  system  makes  Colatz's 
algorithm  useless.  This  is  really  not  surprising  when  one 
considers  that  equation  (91)  is  a  function  of  stochastic 
state  information,  therefore  equation  (91)  is  stochastic. 

To  gain  useful  information  from  this  stochastic  equation,  a 
stochastic  analysis  would  have  to  be  performed.  However, 
note  that  equation  (91)  is  nonlinear  making  its  statistics 
non-Gaussian  and  its  analysis  difficult.  An  alternative  is 
to  use  deterministic  state  information  in  equation  (91)  or 
to  use  another  approach  to  estimating  the  required  step 
size.  Regardless,  this  study  strongly  suggests  that  any 


error  function  used  to  make  decisions  with  regards  to  step 
size  adjustment  must  only  be  based  on  deterministic  infor¬ 


mat ion. 


Tbe  Adams-Moulton  method  is  a  corrector  formula  for 
use  in  a  predictor-corrector  integration  algorithm.  It  is 
probably  the  most  widely  used  corrector  method.  Predictor - 
corrector  formulas  are  iterative  methods  used  to  increase 
the  accuracy  of  solution  at  relatively  low  computational 
costs.  The  Adams-Moulton  formula  requires  derivative  in¬ 
formation  at  the  next  update  increment  in  order  to  compute 
the  state  at  that  same  increment.  The  predictor  formula  is 
used  to  provide  that  derivative  information.  The  predic¬ 
tor-corrector  pair  can  iterate  on  the  solution  for  as  many 
times  as  necessary.  The  iteration  process  can  be  termina¬ 
ted  abruptly  after  some  predetermined  number  of  iterations 
or  it  can  be  terminated  intelligently  after  reaching  some 
defined  convergence  criterion.  The  equation  that  describes 
the  Adams-Moulton  corrector  formula  belongs  to  the  class  of 
multistep  integration  methods.  Thus,  it  only  requires  one 
additional  function  evaluation  per  iteration.  The  specific 
equation  is 

4+1  -  *i  ♦  hC«1Kiiii1>,(i+l)h)  +  ocjKi^ih) 

+  (i-l)h)  +  «4i(x.i_2,<i-2)h)3  (93) 

where  aj  =  9/24;  03  =  19/24;  a3  *  -5/24  ;  0^  =  1/24  ;  and 
k  is  the  number  of  predictor -corrector  iterations.  Note 
that  equation  (93)  is  independent  of  the  method  used  to  ob¬ 
tain  the  prediction  term,  Therefore,  to  analyze 

equation  (93)  on  a  per  iteration  basis,  k  can  be  set  to  1. 


Next,  note  that  the  noise  properties  of  the  corrector  for- 
mula  will  only  be  an  explicit  function  of  four  independent 
noise  terns  derived  at  the  evaluation  of  £(*,*). 

To  analyze  this  system,  the  Jury  method  presented  in 
Chapter  III  will  be  used.  For  this  analysis,  a  difference 
equation  will  be  required.  A  change  of  variable  will  be 
used  to  account  for  the  iterations  in  equation  (93). 

Let  j*i+k*i  +  l.  Using  a  scalar  form  of  equa¬ 
tion  (10)  to  evaluate  f(*,a)  in  equation  (90),  then 
equation  (93)  will  become 

Xj+i  »  Xj_!  +  hCa1(Axj  +  Bwj )  +  cc2(AXj_1  *  BWj.j) 

+  «3(  Axj_2  +  BWj_2^  +  c^(Atj_3  +  (94) 

Taking  the  z-transform  of  equation  (94),  assuming  x0  =  0, 
yields 


zX(z)  *  z-1X(z)  +  hC«x(AX(z)  +  BW(z)) 

+  «2(AX(z)  +  BH( z ) )z“1  +  «3 ( AX ( z )  +  BW(z))z~2 

+  o^(AX(z)  +  BH(  z) )z”33  (95) 

Multiplying  through  by  z3,  collecting  terms,  and  rearrang¬ 
ing  equation  (95),  results  in 

(e*iz3  +  c^z2  +  0C3Z  +  <x4)hBH(z) 

X(z)  *  — - - - - -  (96) 

z4  -  «xhAz3  -  (1  +  012 hA)z2  -  «3hAz  -  o^hA 

Applying  equation  (96)  to  the  Jury  analysis  procedure  given 

in  Chapter  III,  in  the  same  manner  as  equation  (65),  yields 

the  steady-state  expression 


EIx2] 


0.5556(1  -  . 067hA  +  . 012(hA)2)h2B2QD 


(97) 


where 


V(hA) 
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V(hA)  -  -1 . 333hA  -  0 . 6111 ( hA ) 2  +  0.0370(hA)3  +  0.0069(hA)4 

Uncagflljttd  lasflnaigJtsiMa 

V(hA)  in  equation  (97)  should  equal  (1  -  5$).  It  does 
not.  Mhen  equation  (97)  was  first  developed,  the  numerator 
function  was  used  to  determine  the  Qq  to  Q  normalizing  func 
tion  without  considering  the  inconsistency  in  the  denomina¬ 
tor  function.  To  compound  matters,  the  normalizing  func¬ 
tion  worked  well  in  the  simulations.  However,  there  is 
clearly  an  analytical  error  in  this  derivation.  Oblivious 
to  this  error,  the  numerical  analysis  will  be  presented  in 
order  to  show  a  number  of  important  points  regarding  the 
application  of  the  Adams-Noul ton  formula  to  stochastic 
analyses . 

MaaBtlfifll.  JEinfllnaa 

The  numerator  of  equation  (97)  should  equal  B§Qq. 
Assuming  it  does,  equating  it  to  the  scalar  form  of  equa¬ 
tion  (20),  and  evaluating  for  small  hA,  leads  to  the  rela¬ 
tionship 

Qd  =  1 . 8Q/h  (98) 

Figure  17  shows  the  results  of  applying  equation  (98) 
to  a  Adams-Houlton  based  simulation  that  uses  a  4th  order 
Adams-Bashforth  predictor.  To  demonstrate  that  equation 
(98)  is  valid,  independent  of  the  predictor.  Figure  18 
shows  the  results  of  the  Adams-Moulton  simulation  that  uses 
a  4th  order  Runge-Kutta  predictor.  For  the  cases  shown  in 


Figures  17  and  18,  the  algorithm  was  corrected  once  per  up 
date.  Figure  19  shows  the  results  of  the  simulation  using 
the  Adams -Bash forth  predictor  with  4  predictor-corrector 
iterations  per  update. 


Figure  17.  Rdaas-Noulton  Corrector  «/  fldaas-Bashforth  Predictor 
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VI.  SUMMARY  AND  CONCLUSIONS 

This  dissertation  addresses  the  problem  of  determining 
the  correct  relationship  between  the  statistics  of  a  con¬ 
tinuous  random  process  and  the  statistics  of  a  discrete 
random  process  used  to  simulate  the  continuous  random  pro¬ 
cess.  The  findings  of  this  research  are  directly  applica¬ 
ble  to  the  general  field  of  digital  simulation  of  physical 
systems  described  by  ordinary  differential  equations. 

It  is  shown  that  to  ensure  a  faithful  digital  simula¬ 
tion  of  a  continuous  random  process,  the  noise  statistics 
of  the  random  number  generator  must  be  set  to  values  dras¬ 
tically  different  from  the  noise  statistics  of  the  contin¬ 
uous  random  process.  Further,  it  is  established  that  the 
relationship  between  the  continuous  and  discrete  statistics 
will  be  a  function  of  the  integration  method  used  in  the 
digital  simulation. 

The  proper  functional  relationship  between  the  dis¬ 
crete  and  continuous  noise  statistics  was  derived  for 

1.  the  class  of  Runge-Xutta  integrators, 

2.  the  4th  order  Adams -Bash forth  Integrator,  and 

3.  the  Adams-Moulton  corrector  formula. 

Additionally,  the  requirement  for  proper  operation  of  a 
variable-step-size  algorithm  was  developed. 
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The  derived  £unctional  relationship  between  the  dis¬ 
crete  and  continuous  noise  statistics  for  the  class  of 
Runge-Kutta  integrators  results  in  a  unique  normalizing 
factor  for  each  Runge-Kutta  method.  The  derivation  pro¬ 
vides  a  useful  formula  for  calculating  that  factor.  The 
use  of  the  factor  in  the  statistical  relationship  function 
will  provide  the  proper  setting  of  the  random  number  gener 


ator's  statistical  parameters.  The  function  is  demonstrated 
to  be  accurate  for  five  specific  Runge-Kutta  integrators. 

It  is  shown,  in  the  error  analyses  chapter,  that  all 
the  integration  methods  considered  yield  normalization 
functions  that,  without  simplifying  assumptions,  are 
dependent  on  the  system  dynamics.  This  dependency  does, 
however,  decrease  with  the  integration  step  size.  The  rate 
of  decrease  varies  for  each  integration  method.  It  is 
shown  that  the  derived  normalization  functions  for  the 
Runge-Kutta  integrators  are  less  sensitive  to  the  system 
dynamics  than  the  other  examined  methods. 

In  contrast  to  the  Runge-Kutta  methods,  the  derived 
function  for  the  Adams-Bashforth  integrator  is  the  most 
sensitive  to  the  system  dynamics.  The  influence  of  the 
dynamics  can  be  minimized  by  reducing  the  integration  step 
size  to  a  value  much  smaller  than  that  required  for  deter¬ 
ministic  analyses.  However,  the  required  step  size  for 
realizing  an  effective  independence  of  the  system  dynamics 
is  on  the  order  of  step  size  requirements  for  Euler  inte- 


grators  applied  to  deterministic  systems.  This  extremely 
small  step  size  requirement  makes  the  Adams-Bashforth 
integrator  impractical  for  stochastic  simulations. 

The  Adams-Moulton  formula's  normalization  function  is 
derived  independent  of  the  predictor  algorithm.  It  is 
shown  that  the  derived  function  is  valid  when  the  Adams- 
Moulton  corrector  formula  is  used  with  a  Runge-Kutta  pre¬ 
dictor  as  well  as  an  Adams-Bashforth  predictor.  There  is 
an  analytical  inconsistency  in  the  Adams-Mouton  derivation 
that  has  not  been  resolved.  However,  the  numerical  results 
show  a  number  of  useful  findings  for  stochastic  applica¬ 
tions  of  this  often  used  corrector  formula. 

Error  functions  needed  to  implement  variable  step  size 
integration  methods  are  discussed  in  detail.  It  is  shown 
that  when  applied  to  a  stochastic  system,  the  error  func¬ 
tion  for  the  method  considered  will  be  a  stochastic  func¬ 
tion  making  its  usefulness  extremely  limited.  Numerical 
examples  are  presented  to  substantiate  this  finding. 

Though  this  research  investigated  integration  methods 
most  commonly  used  in  practice,  it  is  certainly  not  exhaus¬ 
tive.  However,  the  analysis  procedures  used  in  this  dis¬ 
sertation  are  applicable  to  other  integration  methods  that 
may  be  of  interest.  Without  a  doubt,  the  findings  of  this 
research  establish  the  fact  that  the  proper  relationship 
between  the  statistics  of  the  continuous  random  process  and 
the  statistics  needed  in  the  simulation  for  accurate 


modeling  is  explicitly  dependent  on  the  integration  method 
used  in  the  simulation.  To  ensure  a  faithful  simulation, 
it  is  mandatory  that  the  proper  function  be  derived  and 
validated  for  the  specific  integration  method  used  in  the 
simulation.  Failure  to  do  so  will  likely  invalidate  the 
simulation  results. 
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Application  o£  the  Jury  method  for  analyzing  discrete 
integration  algorithms  of  order  three  and  higher  is  diffi¬ 
cult  without  the  aid  of  computational  tools.  The  difficul¬ 
ty  is  due  to  the  fact  that  the  coefficients  of  the  transfer 
functions  that  describe  the  integration  processes  are  in 
variable  form.  For  instance,  consider  the  transfer  func¬ 
tion  for  the  Adams-Bashforth  integrator  given  by  equation 
(69).  The  coefficients  for  the  numerator  polynomial,  M(z), 
and  the  denominator  polynomial,  L(z),  are 

mg  -  0  lg  «  1 

*1  *  <*ihB  lj  ■  -(1  +  a^hA) 

*2  “  «2hB  I2  “  ~«2hA 

*3  “  «3hB  I3  «  -ajhA 

■4  ■  04 hB  I4  *  -o^hA 

In  the  Jury  algorithm  given  by  equations  (60)  -  (62),  these 
coefficients  will  be  used  to  form  the  matrices  ft  and  ft^  as 
in  equation  (62). 

Recall  from  equation  (61)  that  the  steady-state  value 
of  the  autocorrelation  function  will  be  determined  by 

Thus,  to  solve  the  problem  for  the  4th  order  Adams- 
Bashforth  integrator,  the  determinant  operation  must  be 
accomplished  for  two  fifth  order  square  matrices.  From 
equation  (62),  it  is  seen  that  the  elements  of  ft  will  be  a 
linear  combination  of  the  coefficients  lg  through  I4. 
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Since  the  coefficients  lg  through  I4  are  first-order  poly- 
noaials  in  hA,  the  elenents  of  A  will  also  be  first-order 
polynomials  in  hA.  The  determinant  of  A  will,  therefore, 
be  a  polynomial  in  hA  of  order  five  or  less. 

The  matrix  A^  will  be  the  same  as  A  except  for  the 
first  column.  The  elements  in  the  first  column  of  A}  will 
be  a  combination  of  the  square  of  the  coefficients  mg  -  1114. 
Since  all  the  coefficients  of  M(z)  contain  the  common  term 
hB,  hB  can  be  conveniently  factored  out.  Therfore,  A^  will 
be  a  fifth-order  square  matrix  with  constants  in  the  first 
column  and  first-order  polynomials  in  hA  in  columns  2-5. 
Hence,  the  determinant  of  A}  will  be  a  polynomial  in  hA  of 
order  four  or  less. 

Based  on  the  above  observations,  it  is  expected  that 
the  steady-state  covariance  for  the  Adams -Bash forth  inte¬ 
grator  will  have  the  following  form 

-  ,  „  o  ,  U(hA) 

Etxfl  -  h2B2QD — - - —  «  h2B2QD -  (99) 

lg | A|  V(hA) 

where  U(hA)  is  a  polynomial  in  hA  of  order  four  or  less  and 
V(hA)  is  a  polynomial  in  hA  of  order  five  or  less. 

The  calculation  of  a  determinant  for  a  fifth  order 
matrix  of  constants  is  tedious,  but  for  a  matrix  of  polyno¬ 
mials  it  is  grueling  task.  However,  Jury  provides  a  set  of 
equations  that  perform  the  determinant  operation  for  the 
matrix  given  by  equation  (62).  Though  these  equations  ease 
the  task,  they  only  reduce  the  complexity  slightly.  The 
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equations  are  lengthy  and  notationally  involved.  Because 
of  the  notational  complexity  in  the  equations,  they  were 
re-derived  and  validated  to  ensure  there  were  no  typograph¬ 
ical  errors.  There  were  not.  For  brevity.  Jury's  equa¬ 
tions  will  not  be  reproduced  in  this  document.  However, 
for  those  who  are  really  interested,  the  computer  code 
included  in  this  appendix  contains  a  direct  implementation 
of  Jury's  equations  (see  Listing  1). 

To  determine  U(hA)  and  V(bA),  Jury's  equations  were 
coded  in  a  computer  program  written  in  Pascal.  The  coef¬ 
ficients  for  L(z)  and  M(z)  were  coded  as  a  function  of  the 
program  variable  hA.  The  program  prompts  the  user  for  a 
value  of  hA,  calculates  the  coefficients  of  L(z)  and  M(z), 
evaluates  Jury's  equations,  and  outputs  the  numerical  value 
of  U(hA)  and  V(hA).  Evaluation  of  U(hA)  and  V(hA)  at  six 
distinct  values  of  hA  will  provide  enough  data  to  uniquely 
determine  the  coefficients  of  U(hA)  and  V(hA).  To  provide 
a  little  extra  confidence  that  the  procedure  was  sound, 
seven  points  were  actually  used  knowing  that  if  the  numer¬ 
ical  data  was  valid,  equivalent  results  would  be  obtained. 
For  the  Adams -Bash forth  integrator,  the  results  were 

U(hA)  -  1  -  5 . 6736hA  +  9.5833(hA)2  -  6.25(hA)3  (100) 

V(hA)  -  -2 . 0hA  -  3.333(hA>2  +  7.6736(hA>3 

-9.5833(hA)4  ♦  6.25(hA>5  (101) 

The  question  that  now  arises  is  how  do  these  results 
relate  to  the  function  needed  to  complete  the  analysis; 
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that  is,  BgOp.  Consider,  that  in  steady-state,  the  state 
covariance  equation  is  of  the  form 

ElxflCl  -  3g)  -  BgQD 
which  can  be  rewritten  as 

,  b8od 

Etxf]  *  - —  (102) 

<1  -  3g) 

A  comparison  of  equation  (102)  to  equation  (99)  leads  to 
the  plausible  relationships 

Bg  *  h2B2N(hA)U(hA)  (103) 

and  (1  -  9g)  =  N(hA)V(hA)  (104) 

where  N(hA)  is  a  polynomial  in  hA  that  is  common  to  both  Bg 
and  (1  -  sg).  To  determine  N(hA)  additional  information  is 
required,  such  as  9q. 


The  discrete  state  transition  matrix  maps  the  states 
into  the  states  JLi+i*  Mathematically  this  is  expressed 


li+1  *  5D(h,0)il  (105 

where  Jti  *  jCt^)  and  3Li+j  *  ^.(ti+h).  In  this  derivation, 
all  that  is  known  about  the  specific  problem  will  be  used. 

Consider  the  general  4th  order  Adams -Bash forth  for¬ 
mula,  where 

*4+1  *  Id  ♦  hCoiiKii, tt)  +  «2l<»A-l*ti-l> 

♦  JLi-2' ti-2)  i.i-3*  t^_3)  1  (106 

where  *  55/24;  o«2  ■  -59/24;  03  *  37/24;  and  =  -9/24. 


To  develop  the  specific  problem,  consider 

1.  Restricting  the  problem  to  linear  time-invariant 
differential  systems  given  by 

i<t)  -  A*(t)  +  Bu.(t)  *  f.(jL*t)  (107) 

2.  The  state  transition  matrix  is  independent  of  the 
inputs,  u.(t).  Therefore,  let  jj.(t)  -  0  for  all  t. 

3.  The  continuous  state  transition  matrix  is  given  by 

5(t2,t1)  *  €A<t2  -  tj)  (108) 

The  derivation  of  the  Adams-Bashforth  formula  makes 
the  assumption  that  the  present  and  past  derivatives  are 
known  and  are  exact.  The  same  assumption  will  be  made  in 
this  derivation.  Kith  this  assumption,  equations  (107)  and 
(108)  yield 

I(3Li,ti)  *  A*A 

£(*i-l'ti-l)  *  A*i-1  55  A€_Ah£i 

£<*.i-2'ti-2)  =  A*d-2  =  A€~2Ahii 

i.(*.i-3'ti-3>  *  A*i-3  *  A€_3AhX.i 

Evaluating  equation  (106)  with  these  derivative  functions, 
results  in 

*i+l  *  El  ♦  «iAh  +  «2Ah€“A^ 

+  cc3Ah€-2Ah  +  a4Ah€~3Ah3ii  (109) 

Noting  that  equation  (109)  has  the  same  structure  as  equa¬ 
tion  (105),  it  is  concluded  that 

2q ( h , 0 )  «  I  +  otjAh  ♦  «2Ah€"Ah  ♦  «3Ah€_2Ah  +  c^AhC”3^  (110) 

Notice  that  unlike  the  Runge-Kutta  methods,  the  dis¬ 
crete  state  transition  matrix  for  the  Adams-Bashforth  inte- 
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grator  is  an  infinite  series.  The  question  arises,  does 
this  state  transition  Matrix  approximate  the  Taylor  series 
fora  of  3(h,0)?  To  answer  this  question,  expand  equation 
(110)  in  a  series  fora  , collect  teras  and  coapare  it  to 
3(h,0)  -  €hA 

-  I  ♦  Ah  ♦  (Ah)2/ 2  ♦  ( Ah)3/6  +  (Ab)4/24 

♦  (Ah)*/120  +  ...  (Ill) 

Expanding  the  exponentials  in  equation  (110)  provides 
€“hA  «  I  -  hA  +  ( h A ) 2 / 2  -  <hA)3/6  ♦  (hA)4/24  ■»■... 

€“2hA  «  i  -  2hA  +  2(hA)2  -  4(hA)3/3  ♦  2(hA)4/3  ♦  ... 
€-3bA  »  I  -  3hA  ♦  9(hA)2/2  -  9(hA)3/2  ♦  27(hA)4/8  +  ... 
Evaluating  equation  (110)  using  these  expansions  and  the 
oil's  given  in  equation  (106),  results  in 

2D(h,0)  »  I  ♦  Ah  ♦  (Ah)2/ 2  +  (Ah)3/6  +  (Ah)4/24 

-  49(Ah)5/120  ♦  ...  (112) 

A  coaparison  of  equation  (112)  to  equation  (111)  shows  that 
3q(A,0)  s  3(h,0)  with  error  on  the  order  of  h5.  With  3q  in 
hand,  the  Jury  analysis  can  be  completed. 

Determination  of  H(hA)  and  Bg 

The  general  form  of  N(hA)  is  assumed  to  be 

N(hA)  -  *g  +  YjhA  +  Y2 (hA ) 2  +  Y3 ( hA ) 3  ...  (113) 

The  coefficients  in  equation  (113)  can  be  computed  di¬ 
rectly  from  equation  (104).  Consider  that 

1  -  2§  -  -2hA  -  2 ( hA ) 2  -  4(hA)3/3  -  2(hA)4/3  ♦  ...  (114) 

Talcing  the  product  of  equation  (113)  with  equation  (101) 
and  equating  the  resulting  coefficients  to  the  coefficients 
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in  equation  (114)  results  in  the  following  constraint  equa¬ 
tions 

-2*g  -  -2 

-3. 3333*g  -  2*x  -  -2 
7 . 6736  Vg  -  3.3333*!  -  2 *2  11  -4/3 
-9.5833*g  +  7.6736*!  -  3.3333*2  -  2*3  -  -2/3 
and  so  on.  Solving  for  *g  through  *3  results  in 

N(hA)  -  1  -  2<hA)/3  +  5.61(hA)2  -  16.379(hA)3  +  ...  (115) 


Solving  for  Bg  in  equation  (103)  results  in 

Bg  -  h2B2(l  -  6. 34hA  +  18. 98(hA)2  -  60.84(hA)3  +...)  (116) 


',rVn.r\  rtix  * 


Listing  1. 

PASCAL  Prograa  Used  to  Evaluate  Jury's  Equations 


Prograe  Jury_Deterainant; 

(  Calculates  U(hA)  and  V(hA)  for  the  Adaas-Bash forth 
integrator.  The  coefficients  10  through  14  have  been 
replaced  by  a0  through  a4  to  enhance  readability  ) 

var 

hA,a0,al,a2,a3,a4,el,e2,e3,e4,e5, 
alfal,alfa2,al fa3,alfa4, 
a0,al,a2,a3,a4,B0,Bl,B2,B3,B4  :  real; 

Q0»Q1*Q2/Q3,Q4,BQ, AQ  :  real; 

begin 

alfal  s-  55.0/24.0; 
alfa2  -59.0/24.0; 

alfa3  37.0/24.0; 

alfa4  :■  -9.0/24.0; 

(  Proapt  User  for  input  and  Get  It  } 
write ('Input  hA  ');  readln(hA); 

(  Define  Coefficients  ) 

a0  ;■  0.0; 
al  :■  alfal; 
a2  :■  alfa2; 
a3  :■  alfa3; 
a4  : ■  alfa4; 
a0  i . 0 • 

al  -(1^0  4  alfal*hA); 

a2  :«  -alfa2*hA; 
a3  : ■  -alfa3*hA; 
a4  :■  -alfa4*hA; 

<  Jury's  Equations  } 

el  :■  a0  4-  a2; 
e2  :■  al  4  a3; 
e3  :>  a2  4-  a4; 

e4  :•  a0  4-  a4; 
e5  : ■  a0  4  a2  4  a4; 

B0  :■  a0*a0  4  al*al  4-  a2*a2  4  m3*a3  4-  a4*a4; 

B1  :«  2.0*(a0*al  4  al*a2  4  a2*a3  4  a3*a4); 

B2  2.0*(a0*a2  4  ml*a3  4  a2*a4); 

B3  :•  2.0*(a0*a3  4  al*a4); 

B4  s-  2.0*a0*a4; 

00  :■  a0*(el*e4  -  a3*e2)  4  a4*(al*e2  -  e3*e4); 

01  a0*(al*e4  -  a2*a3)  4  a4*(al*a2  -  a3*e4); 

02  : ■  a0*(al*e2  -  a2*el)  4  a4*(a2*e3  -  a3*e2); 

03  al*(al*e2  -  e3*e4)  -  a2*(al*el  -  a3*e3> 

♦  —  a3*a9)i 

04  ;»  a0*(e2*(al*a4  -  a0*a3)  4  e5*(a0*a0  -  a4*a4)) 

4  (e2*e2  -  e5*e5)*(al*al  -  al*a3  4  (a0  -  a4)*(e4  -  a2)) 


Continuation  of  Listing  1. 

{  BQ  is  U(hA)  and  AQ  is  V(hA)  ) 

BQ  :>  a0*B0*Q0  -  a0*Bl*Ql  +  a0*B2*Q2  -  a0*B3*Q3  +  B4*Q4; 
AQ  a0*( (a0*a0-a4*a4) *Q0  -  (a0*al  -  a3*a4)*Ql 
+  (a0*a2  -  a2*a4)*Q2  -  (a0*a3  -  al*a4)*Q3); 

(  Output  to  terainal  } 
writeln( *Por  hA  *  ',hA); 
writeln('Q0  -  ',Q0, '  Q1  -  ',Q1>; 
writeln('Q2  *  *,02,*  Q3  -  »,Q3,'  Q4  -  »,Q4); 
writelnt'AQ  -  *  ,AQ,  '  BQ  -  ' , BQ ) ; 

(  Output  to  Printer  } 
writeln( 1st); 

writelndst, 'For  hA  *  ',hA); 
writelndst,  'Q0  -  ',Q0,  »  Q1  -  »,Q1); 

writelndst,  '02  =  ',Q2, '  Q3  **  *,Q3,  '  Q4  *  ',Q4); 

writelndst,  *AQ  «  ' ,  AQ, '  BQ  *  ' , BQ ) ; 

writelndst); 
writelndst); 
end. 
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APPENDIX  B 


CODE  FOR  RANDOM  NUMBER  GENERATOR 


USED  IN  SIMULATION 


u  u  u 


»*« 


100 

Listing  2. 

FORTRAN  Code  o£  Random  Number  Generator 
with  Gaussian  Shaping  Algorithm 


REAL  FUNCTION  GAUSSN(SIG,SEED)  I 

« 

**  FOR  GOOD  RESULTS  USE  AN  INITIAL  SEED  *  824064364  j 

INTEGER *4  SEED 

GNOIZ*0.  j 

DO  10  1*1, 12  ! 

GNOIZ  *  GNOIZ  +  URAND( SEED) 

10  CONTINUE  t 

GAUSSN«SIG*( GNOIZ -6.0)  \ 

RETURN  t 

END  fl 

REAL  FUNCTION  URAND(SEED)  { 

INTEGER*4  B2E15,B2E16, MODLUS, HIGH15,HIGH31,LOH15, LOHPRD, 

6  MULTI , MULT2 , 0 VFLOH , SEED 
DATA  MULTI, MULT2/24112, 26143/ 

DATA  B2E15,B2E16,MODLUS/32768, 65536, 2147483647/  | 

HIGH15  *  SEED/B2E16  •' 

LOHPRD  *  (SEED  -  HIGH15*B2E16 ) *MULT1  i 

LOH1S  *  LOHPRD/B2E16 

HIGH31  *  HIGH15*MULT1  +  LOW 15  l 

OVFLOH  *  HIGH31/B2E15  I 

SEED  *  (  (  ( LOHPRD  -  LOW15*B2E16)  -  MODLUS)  +  !• 

6  (HIGH31  -  OVFLOW*B2E15 )  *B2E16 )  ♦  OVFLOH  \\ 

IF  ( SEED .  LT .  0 )  SEED  -  SEED  +  MODLUS  L; 

HIGH15  *  SEED/B2E16  j 

LOHPRD  *  (SEED  -  HIGH15*B2E16 ) *MULT2  ! 

LOH15  -  LOHPRD/B2E16  * 

HIGH31  *  HIGH15*MULT2  +  LOH15  l 

OVFLOH  -  HIGH31/B2E15  K 


i 


SEED  *  (((LOHPRD  -  LOH15*B2E16)  -  MODLUS)  + 

6  (HIGH31  -  OVFLOH*B2E15)*B2E16)  +  OVFLOH 

IF  ( SEED . LT . 0 )  SEED  *  SEED  +  MODLUS 

URAND  -  FLOAT ( 2* (SEED/ 256)  +  1 )/ 16777216 . 0 

RETURN 

END 


